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ABSTRACT 

We derive conservative, multidimensional, energy-dependent moment equations for neutrino transport in 
core-collapse supernovae and related astrophysical systems, with particular attention to the consistency of con- 
servative four-momentum and lepton number transport equations. After taking angular moments of conserva- 
tive formulations of the general relativistic Boltzmann equation, we specialize to a conformally flat spacetime, 
which also serves as the basis for four further limits. Two of these — the multidimensional special relativistic 
case, and a conformally flat formulation of the spherically symmetric general relativistic case — are given in ap- 
pendices for the sake of comparison with extant literature. The third limit is a weak-field, 'pseudo-Newtonian' 
approach (Kim et al. 2009, 2012) in which the source of the gravitational potential includes the trace of the 
stress-energy tensor (rather than just the mass density), and all orders in fluid velocity v are retained. Our pri- 
mary interest here is in the fourth limit; '©(w)' moment equations for use in conjunction with Newtonian self- 
gravitating hydrodynamics. We show that the concept of '0{vy transport requires care when dealing with both 
conservative four-momentum and conservative lepton number transport, and present two self-consistent op- 
tions: 'C'(i')-plus' transport, in which an 0{v^) energy equation combines with an 0{v) momentum equation 
to give an 0{v'^) number equation; and '©(w) -minus' transport, in which an 0{v) energy equation combines 
with an 0(1) momentum equation to give an 0{v) number equation. 

Subject headings: neutrinos — radiative transfer — supernovae: general 



1. INTRODUCTION 

Detailed simulations of neutrino transport are central to 
understanding the core-collapse supernova explosion mech- 
anism (e.g., Mezzacappa 2005; Woosley & Janka 2005; 
Kotake et al 2006; Janka 2012). Neutrinos of all flavors (i.e., 
electron, mu, and tau neutrinos and antineutrinos) carry away 
~ 99% of the gravitational energy released (~ 10^'^ erg) dur- 
ing core collapse of a massive star (M > 8 Mq). As the 
dominant energy source (in the case of slowly-rotating pro- 
genitors), it is likely that the neutrinos streaming from the 
neutrinospheres of a proto-neutron star (PNS; the collapsed 
core of a massive star) deposit enough energy into the stellar 
fluid to power the supernova: about5-10% (~^ 10^^ erg) of the 
neutrino energy emitted in the first second after core bounce 
(> 10^2 erg) must be transferred to give rise to a neutrino- 
driven supernova. 

Certainly, robust and accurate numerical methods must be 
deployed. Numerical simulations involving multi-physics 
codes (coupling neutrino radiation transport, magnetohydro- 
dynamics, nuclear reaction kinetics and equations of state, 
and gravity, preferably — and ultimately — in general relativ- 
ity) must be evolved for ^ 0(10^-10^) timesteps to study in 
detail the energy transfer from the neutrino radiation field to 
the supernova matter, and the onset and development of the 
explosion. In particular, conserved quantities (i.e., total en- 
ergy and lepton number) must be preserved within tolerable 
margins, before firm conclusions can be drawn from the sim- 
ulation outcomes. 

So far, only spherically symmetric simulations have 
reached the desired level of realism in neutrino transport, 
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where general relativistic, angle- and energy-dependent neu- 
trino radiation hydrodynamics simulations can be performed 
(e.g., Liebendorferet al. 2001, 2004; Sumiyoshi et al. 
2005). However, spherically symmetric simulations do not 
result in neutrino-driven explosions (Rampp & Janka 2000; 
Liebendorfer et al. 2001; Thompson et al. 2003, except in 
the case of the lightest progenitors with O-Ne-Mg cores; 
Kitaura et al. 2006), and, based on results from 2D (ax- 
isymmetric) simulations, there is now an emerging consen- 
sus that multidimensional effects due to hydrodynamics insta- 
bilities are important in shaping the core-collapse supernova 
explosion (e.g., Herantetal. 1994; Burrows et al. 1995; 
Janka & Muller 1996; Fryer & Warren 2002; Burrows et al. 
2006, 2007; Bruenn et al. 2009; Marek& Janka 2009; 
Suwaetal. 2010; Muller etal. 2012a,b; Bruenn etal. 
2012). 

Neutrino transport in multidimensional simulations is far 
less mature. Discrepant results from 2D simulations by dif- 
ferent research groups, in part due to the physics included 
in the models and the approximations made, remain to be 
fully investigated. It appears clear, however, that energy- 
dependent neutrino transport (including "observer correc- 
tions" due to fluid motions and inelastic scattering; e.g., 
Lentz et al. 2012a) as well as a comprehensive set of weak- 
interactions are required (e.g., Lentz et al. 2012b). 3D sim- 
ulations are still in their infancy, although recent progress 
has been made (e.g.. Fryer & Young 2007; Takiwaki et al. 
2012). 

Current approaches to neutrino transport in multidimen- 
sional simulations include the so-called ray-by-ray ap- 
proximation with the multigroup (energy-dependent) flux- 
limited diffusion approximation (e.g., Bruenn et al. 2009, 
2012), the two-moment model with approximate Boltz- 
mann closure (e.g., Marek & Janka 2009; Muller etal. 
2012a), and the so-called isotropic diffusion source approx- 
imation (e.g., Liebendorfer et al. 2009; Suwaetal. 2010). 
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Other recent approaches to neutrino radiation transport in- 
clude multidimensional, multigroup flux-limited diffusion — 
vaUd to 0{v) (e.g., Swesty & Myra 2009; Zhang et al. 
2012) or to 0(1) (e.g.. Burrows et al. 2006)— or multi- 
dimensional, energy-integrated (grey) or multigroup two- 
moment models with analytic closure (e.g., Kuroda et al. 
2012; Obergaulinger & Janka 2011; Shibata et al. 2011). 
Ott et al. (2008) simulated the post-bounce phase of core- 
collapse supernovae in axial symmetry using 0(1) multi- 
group and multiangle neutrino transport (see also the recent 
developments by Sumiyoshi & Yamada 2012). 

Lack of sufficient computational resources is the primary 
reason for the various approximations made in multidimen- 
sional simulations of neutrino transport in core-collapse su- 
pernovae. Gravitational collapse, a steep density cliff at 
the contracting surface of the PNS, and regions of post- 
shock turbulence recommend regions of high spatial res- 
olution. Moreover, differences in the evolution of im- 
portant hydrodynamic phenomena in two and three spatial 
dimensions — for example, the Standing Accretion Shock In- 
stability (e.g., Blondin et al. 2003; Blondin & Mezzacappa 
2007), or turbulence (Ishihara et al. 2009; Boffetta & Ecke 
2012) — suggest that simulations should ideally (and eventu- 
ally) be performed in 3D. Recent explorations using so-called 
neutrino light-bulb models emphasize differences between 2D 
and 3D simulations (e.g., Nordhaus et al. 2010; Hanke et al. 
2012; Burrows et al. 2012; Couch 2012). Models based on 
evolving angular moments of the neutrino distribution func- 
tion are therefore attractive for 3D simulations, where the so- 
lution space must be reduced. 

With an eye toward future multidimensional simula- 
tions of angle- and energy-dependent neutrino transport, 
Cardall & Mezzacappa (2003, hereafter CM03) derived con- 
servative general relativistic kinetic equations intended for 
implementation in numerical codes for simulations of neu- 
trino transport in core-collapse supernovae. They adopted dis- 
tinct position and momentum space coordinates, by employ- 
ing spacetime position coordinates associated with a global 
coordinate basis and momentum space coordinates associated 
with an orthonormal comoving basis. The kinetic equations 
derived by CM03 describe the evolution of the specific par- 
ticle number density and four-momentum (six-dimensional 
phase space densities). Such conservative equations allow 
for convenient treatment of neutrino-matter interactions in the 
comoving-frame, and may be particularly useful when con- 
structing numerical methods for the multidimensional Boltz- 
mann equation with desirable conservation properties (i.e., 
conserves total lepton number and energy). However, Boltz- 
mann neutrino transport in 3D core-collapse supernova simu- 
lations will not be possible for some time to come (exascale 
computational resources will be required, even with moderate 
phase space resolution), and approximate methods will con- 
tinue to play a dominant role in the foreseeable future. 

In preparation for development of numerical methods for 
neutrino radiation transport in our astrophysical simulation 
code GenASiS (Cardall et al. 2012a,b), we derive conserva- 
tive, multidimensional, monochromatic moment equations for 
neutrino radiation transport from first principles, employing 
the general relativistic framework of CM03. We derive a two- 
moment model, which evolves angular moments of the neu- 
trino distribution function; i.e., the neutrino energy and mo- 
mentum (as opposed to only the energy in the flux-limited 
diffusion approach). The radiation moments are functions 
of spacetime position components associated with a global 



coordinate basis x^^ and the radiation energy measured by a 
comoving observer e. First we specialize the moment equa- 
tions to the case of a conformally flat spacetime. Then we 
specialize the equations further to the pseudo-Newtonian (cf. 
Kim et al. 2012) and the Newtonian gravity, 0{v) limits. 
Special relativistic equations are presented in Appendix A, 
while we present general relativistic moment equations for 
spherically symmetric spacetimes in Appendix B. (General 
relativistic moment equations employing the 3h-1 formula- 
tion have been presented separately in Cardall et al. 2012; 
see also Shibata et al. 2011.) We elucidate the relationship 
between the moment equations for neutrino four-momentum 
and neutrino number in general (see also Cardall et al. 20 1 2), 
and we pay special attention to consistency between the con- 
servative moment equations for neutrino four-momentum and 
the conservative equation for the neutrino number density in 
the non-relativistic limits. Detailed knowledge of this rela- 
tionship may be useful when constructing consistent numer- 
ical methods for neutrino radiation transport, where lepton 
number conservation (in addition to total energy conserva- 
tion) must also be considered. Our conservative moment 
equations correspond to similar non-conservative moment 
equations presented by other authors (e.g., Buchler 1979, 
1983, 1986; Kaneko et al. 1984; Munier & Weaver 1986b; 
Mulleretal. 2010). 

The 0{v) limit of the radiation moment equations is not 
uniquely defined, and care should be taken to ensure con- 
sistency between the conservative equations for the radiation 
four-momentum and the neutrino number equation (see also 
Cardall et al. 2005, for a detailed discussion). In the non- 
relativistic limit, we find that exact consistency between the 
conservative moment equations for neutrino four-momentum 
and the conservative neutrino number equation can be ob- 
tained by adopting different orders of v for the radiation en- 
ergy and momentum equations, respectively. As an illustra- 
tive example, in the laboratory frame the neutrino number 
density £j^ and the neutrino energy density £ and momentum 
density J^' are to 0{v) related by 

e£j^^£-v,J^\ (1) 

where Vi is the fluid three-velocity. The evolution equations 
for the respective quantities are similarly related. By approxi- 
mating the conservative (lab-frame) radiation energy and mo- 
mentum equations to 0{v) and 0(1), respectively, we ob- 
tain exact consistency with the conservative 0{v) neutrino 
number equation, and we avoid the introduction of 'extra' 
O(w^) terms. However, because we use an 0(1) momentum 
equation, this description is formally only accurate to 0(1) 
(we call this the 0(w)-minus approximation to indicate that 
velocity-dependent terms are included, and to distinguish it 
from 'traditional' 0(1) descriptions, where no moving fluid 
effects are included). We obtain a description that is formally 
accurate to 0{v) by adopting 0{v^) and 0{v) radiation en- 
ergy and momentum equations, respectively. This description 
is exactly consistent with the conservative 0{v'^) neuterino 
number equation. We refer to this system of moment equa- 
tions as the 0(u)-plus radiation moment equations. 

We have organized the paper as follows: we briefly summa- 
rize general relativistic kinetic equations in Section 2. Angu- 
lar moments of the distribution are defined in Section 3, and 
we derive conservative evolution equations for the so-called 
zeroth- and first-order angular moments in Section 4. (We de- 
rive conservative monochromatic equations for the neutrino 
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number density and four-momentum.) We also elucidate the 
relationship between the stress-energy equation and the num- 
ber equation in Section 4. In Section 5 we specialize the 
moment equations for a sufficiently general (conformally flat) 
spacetime metric to accommodate the pseudo-Newtonian ra- 
diation moment equations; the Newtonian gravity, 0{v) equa- 
tions; the special relativistic moment equations; and the gen- 
eral relativistic moment equations for spherically symmetric 
spacetimes. The pseudo-Newtonian moment equations are 
presented in Section 6. We introduce the Newtonian grav- 
ity, 0(u)-plus and 0(ij)-minus approximations in Section 7. 
A summary and discussion is given in Section 8. Special rela- 
tivistic moment equations are presented in Appendix A. Gen- 
eral relativistic moment equations for spherically symmetric 
spacetimes are presented in Appendix B, where we also com- 
pare the equations with those solved by Miiller et al. (2010). 

2. GENERAL RELATIVISTIC KINETIC EQUATIONS 

In this section we briefly summarize general relativistic ki- 
netic equations, and define the variables associated with the 
equations (see for example Lindquist 1966; Ehlers 1971; 
Israel 1972, for detailed treatments of relativistic kinetic the- 
ory). The kinetic equations form the basis for our derivation 
of angular moment equations for neutrino radiation transport. 
To this end, we consider a single species of massless, elec- 
trically neutral particles (i.e., classical neutrinos). We adopt a 
'geometrized' unit system in which the vacuum speed of light, 
the gravitational constant, and the Planck constant are unity, 
and we adopt the usual Einstein summation convention and 
let repeated Greek indices run from to 3, and repeated Latin 
indices run from 1 to 3. 

The classical neutrino distribution function f{x^,p^) — the 
phase space particle density — is governed by the Boltzmann 
equation. In its geometric form, it states that the change in 
/ along a phase space trajectory, parametrized by A, equals 
the density C [/] of point-like collisions that add or remove 
particles from the trajectory; i.e.. 



dX 



<C[f] 



(2) 



For practical computations it is necessary to introduce phase 
space coordinates. By introducing spacetime coordinates x^, 
and momentum space coordinates p\ Equation (2) becomes 



dxf" df dp^ df 



= C [/] . 



(3) 



(By the mass shell constraint p/i = 0, only three of the 
four-momentum components are independent, and the con- 
traction over momentum coordinates in Equation (3) only 
runs over three indices.) The geodesic equations describing 
the trajectory can be written as 

(4) 



d\ 

d\ ^ 



(5) 



Then, the general relativistic Boltzmann equation reads (e.g., 
Lindquist (1966); Mezzacappa & Matzner (1989); CM03) 



up 



9p* 



'[/] 



(6) 



In the above we have used our freedom in choosing dis- 
tinct spacetime and momentum space coordinates; {a;''} = 



{t, a;'} are spacetime position components in a global coor- 
dinate (holonomic) basis, while {p^} = e{l, n*} are four- 
momentum components in a comoving orthonormal basis, 
where n' are components of the spatial unit vector parallel 
to p', with basis tangent to the spatial coordinate lines in a 
comoving orthonormal basis. 

A 'comoving-frame' at a spacetime point is an inertial ref- 
erence frame whose four-velocity instantaneously coincides 
with that of a moving fluid element at that point, and is the 
frame where interactions between the radiation field and the 
fluid are most easily described (Mihalas & Mihalas 1999). 
In the comoving frame, the four-velocity of a comoving ob- 
server is {u^} ~ {1,0,0,0}. Quantities with unaccented 
indices are defined with respect to the global coordinate ba- 
sis, while quantities whose indices are accented with a hat 
(comoving-frame quantities) are defined with respect to the 
orthonormal comoving basis. The composite transformation 
L'i = e'^-j A^. transforms four-vectors defined with respect 
to the orthonormal comoving basis into four-vectors defined 
with respect to the global coordinate basis. Thus, the four- 
velocity of a comoving observer with respect to the global 
coordinate basis is = L''^ u^-. The composite transfor- 
mation consists of a Lorentz transformation A^- from the or- 
thonormal comoving basis (comoving-frame) to an in general 
noncomoving orthonormal tetrad basis (the orthonormal lab- 
frame; e.g., ~ A^^ u^), followed by a transformation e^^ 
from the orthonormal tetrad basis to the global coordinate 



basis (e.g., 



e n u^). Quantities with indices accented 



with a bar (lab-frame quantities) are defined with respect to 
the orthonormal tetrad basis. The "tetrad transformation" lo- 
cally transforms the metric tensor g^,y into the Minkowski 
metric tensor; e'^^e^'pg^^ = rjf^i^ = diag[-l, 1, 1, 1]. (The 
Minkowski metric is preserved by the Lorentz transforma- 
tion, and the composite transformation is also a tetrad trans- 
formation.) The inverse composite transformation is denoted 
L'), = A*^^ e'^; i.e., L'i L^^ ^ 5^'^, where is the Kronecker 
delta. 

The 'connection coefficients' associated with the orthonor- 
mal comoving basis T^.. (cf. Equation (6)) are expressed in 
terms of the connection coefficients associated with the global 
spacetime coordinate basis as 



^M^P Q^p 



(7) 



where the coordinate basis connection coefficients are given 
in terms of the spacetime metric g^i^ as 



1 



per 



^ dg^p _ dg^ 
\ dxP dx"" dx° 



(8) 



and g^'" is the contravariant metric tensor; g^'' = 5^-^. 
Note now that the factors outside the momentum space deriva- 
tive of / in the second term on the left-hand side of Equation 
(6) are proportional to space and time derivatives of the fluid 
three-velocity and components of the metric tensor These 
terms account for Doppler and Einstein shifts. See for ex- 
ample Lentz et al. (2012a) for a general discussion of these 
terms and a demonstration of their importance in simulations 
of neutrino transport in core-collapse supemovae. 

The collision operator C [/] describes the momentum space 
evolution of the distribution function / due to neutrino-matter 
interactions. Neutrino-matter interactions crucial to realistic 
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simulations of CCSNe include particle creation and destruc- 
tion via emission and absorption, scattering (both elastic and 
inelastic), and neutrino pair creation and annihilation. We will 
not discuss neutrino-matter interactions in any detail in this 
paper See, however, Baienn (1985); Bun^ows & Thompson 
(2004), and the references therein for general details. See also 
Lentz et al. (2012b) for a recent discussion of neutrino opac- 
ities in core-collapse supernova simulations. 

Our choice of phase space coordinates is well suited for 
numerical solution of the Boltzmann equation (or equations 
based on angular moments of the Boltzmann equation). In 
Equation (6), the particle distribution function is parametrized 
in terms of spacetime position components in a global coordi- 
nate basis {x^} and momentum components in an orthonor- 
mal basis comoving with the fluid {p^}. This specific choice 
is motivated by our intent to develop numerical methods for 
computer simulations of neutrino transport. Neutrino-matter 
interactions are on the one hand most easily handled computa- 
tionally in the frame comoving with the fluid, where material 
properties are isotropic (e.g., Mihalas & Mihalas 1999). On 
the other hand, when integrated over momentum space, the 
Boltzmann equation expresses conservation of particle num- 
ber and energy in the laboratory frame. Accurate account- 
ing of total lepton number and energy in numerical simula- 
tions of core-collapse supemovae is important and extremely 
challenging. Convenient treatment of neutrino-matter interac- 
tions and global conservation are naturally expressed in this 
phase space coordinate basis. However, Equation (6) is not 
in conservative form, and, from a practical standpoint, may 
not be the best starting point for designing numerical meth- 
ods to be implemented in simulation codes. Numerical meth- 
ods based on the solution of the Boltzmann equation as it is 
expressed in Equation (6) will in general not conserve lepton 
number or energy. Even for numerical methods based on a 
number-conservative formulation of the Boltzmann equation, 
care must be taken in the discretization process in order to 
ensure total lepton number and energy conservation within 
acceptable limits, so that firm conclusions can be drawn from 
the simulations (Liebendorfer et al. 2004). 

To help facilitate the realization of lepton number and en- 
ergy conservation in multidimensional simulations of neu- 
trino transport, CM03 derived conservative general relativis- 
tic formulations of kinetic theory. In particular, their number- 
conservative reformulation of Equation (6) (cf. Eq. (168) in 
CM03) reads 



where the spacetime divergence is 

1 d 



and the momentum space divergence is 

dp 



(9) 



(10) 



_d_ 



(1- 



IP 

det 



dct 
dp 



du 



du 



(11) 



The determinant of the metric tensor is denoted g. 
Similarly, the four-momentum-conservative reformulation of 
Equation (6) (cf. Eq. (169) in CM03) reads 



§^[/]+M^[/] = L"-/C[/], 



(12) 



where the spacetime divergence is 
1 d 



' —g dx 
+T^^,L%L^,pPp^ f 

and the momentum space divergence is 

dp^ ^1 



(13) 



M^[/] = -|p 
d ( 1 



det 



det 
dp 



du 



du 



du 



r^.g^L'^.p'p'p'f)- (14) 



Equations (9) and (12) are conservative in the sense that, 
when integrated over momentum space, the momentum space 
derivatives vanish, and the resulting equations — expressing 
number and four-momentum balance, respectively — are fa- 
miliar from position space conservation law theory (see 
CM03, for a detailed discussion). Note that Equations (9) 
and (12) are merely conservative (nontrivial) reformulations 
of Equation (6). Thus, Equations (9) and (12) are, of course, 
not independent. While a numerical method is based on solv- 
ing one of the equations, the numerical solution should ideally 
remain consistent with both. A formidable challenge for fu- 
ture efforts will be to construct a discrete representation of the 
multidimensional Boltzmann equation (Equation (9) or Equa- 
tion (12)) that is both number and energy conservative (i.e., si- 
multaneously consistent with Eqs. (9) and (12)) — as was done 
by Liebendorfer et al. (2004) in the spherically symmetric 
case. The conservative formulations provided by Equations 
(9) and (12) may be helpful for this task. Numerical meth- 
ods based on solving for angular moments of the distribution 
function should also retain this "dual consistency." In this pa- 
per we will discuss requirements and prospects for construct- 
ing numerical methods for neutrino transport in CCSNe based 
on moment models that simultaneously conserve total lepton 
number and total energy. 

In Equations (9) and (12), a change to spherical momen- 
tum space coordinates {u^} = {e, ip} has been facili- 
tated. The Cartesian momentum components are expressed 
in terms of the spherical momentum space coordinates by 

{p^} = i{n^, rP' , TV'} — e{cost?, sint?cos(y3, sin Tasini/?}. 
The Jacobian matrix associated with the transformation is 



dp^ 



cos-i? 



-esint? 







— — : = sin?9cosiy9 e CDS'!? cos — esirn?sin(^ (15) 



du^ 



sin I? sin ip e cos •& sin ip e sin d cos ip 



whose inverse is 

du' _ 1 
dpi e 




e sin "d cos ip e sin -d sin p> \ 
cos 1? cos (/3 cos sin I . (16) 
— svtvip/ sin •§ cos ip/ sin / 



The particle energy measured by a comoving observer (as- 
suming massless neutrinos) is |p| = e, and the determinant of 
the Jacobian matrix (dp^/du^) is 



det 



dp 



du 



■ sin??. 



(17) 



This concludes our summary of general relativistic kinetic 
equations. 

3. ANGULAR MOMENTS OF THE DISTRIBUTION FUNCTION 
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The distribution function /(a-^,p*) provides a detailed sta- 
tistical description of the particle momentum at a given space- 
time location. However, as a seven-dimensional object it 
is in general computationally prohibitive to compute with 
sufficient phase space resohition in multidimensional super- 
nova simulations using currently available computer hard- 
ware. The solution space must be reduced. Instead of solv- 
ing the Boltzmann equation directly for the particle distri- 
bution function, solving equations for moments of the dis- 
tribution function has become a popular method of reducing 
the dimensionality of the transport problem to a computation- 
ally tractable one (e.g., Lindquist 1966; Anderson & Spiegel 
1972; Castor 1972; Thorne 1981; Munier & Weaver 1986b; 
Mihalas & Mihalas 1999). Variants of so-called moment 
models (with various degrees of sophistication) have been 
used extensively for simulation of neutrino transport in core- 
collapse supernovae (e.g., Bruenn 1985; Herant et al. 1994; 
Rampp & Janka 2002; Bruenn et al. 2006; Buras et al. 
2006; Burrows etal. 2006, 2007; Bruenn et al. 2009; 
Marek & Janka 2009; Swesty & Myra 2009; Miiller et al. 
2010; Obergaulinger & Janka 2011; Bruenn etal. 2012). 
Energy-dependent, general relativistic angular moment equa- 
tions are presented in this paper (see also Miiller et al. 2010; 
Shibataetal. 201 1; Cardall et al. 2012). 

The neutrino heating rate in the gain region (i.e., the re- 
gion behind the shock where the neutrino heating rate exceeds 
the neutrino cooling rate) depends quadratically on the neu- 
trino energy spectrum (e.g., Mezzacappa 2005), and mainly 
for this reason do we retain the energy dependence, while we 
integrate over the momentum directions. The dimensionality 
of the transport problem is then reduced to five (four space- 
time dimensions, and one momentum space dimension). In 
particular, we wish to derive equations for the lab-frame neu- 
trino number density and four-momentum density expressed 
in terms of angular moments of the distribution function de- 
fined in the comoving-frame. The invariant momentum space 
volume element is 



dr'u 



esind deM dip. 



(18) 



It is expressed in terms of spherical comoving-frame momen- 
tum space coordinates on the far right side. 

The number-flux four- vector (integrated over the entire mo- 
mentum space Vp) is defined in terms of moments of the dis- 
tribution function as (e.g., Lindquist 1966) 



d^p 



de. 



(19) 



where the comoving-frame neutrino number density and num- 
ber flux density are denoted Jj^ and T-L ^, respectively. We use 
the calligraphic font to distinguish energy-dependent radia- 
tion quantities from energy-integrated (grey) radiation quan- 
tities (e.g., vs. A^^). Note that what we here call the 
lab-frame number-flux four-vector is a function of the co- 
ordinate basis spacetime position components and the 
neutrino energy measured by a comoving observer e; i.e., 
~ A/'^(a;^,e). This parametrization of radiation trans- 
port variables was suggested by CM03 (see also Riffert 1986; 
Mezzacappa & Matzner 1989). This is also the parametriza- 
tion used by Shibata et al. (201 1) and Cardall et al. (2012). 

In a similar manner, the momentum space integrated sym- 
metric stress-energy tensor is defined in terms of moments of 
the distribution function as (e.g., Lindquist 1966) 

T^^^ [ p^^p-fpE^^ r T^- de, (22) 

Jv^ \P\ ^ Jo 

where the monochromatic stress-energy tensor in the 
comoving-frame is defined in terms of angular moments of 
the distribution function as 



1 



pf'p" fdn. 



(23) 



The monochromatic lab-frame stress-energy tensor T^" = 
L^^ L^P T^'^ is also a function of the neutrino energy mea- 
sured by a comoving observer; i.e., T'^" = T^'^{x^ , e). Fur- 
thermore, we define the components of the monochromatic 
comoving-frame stress-energy tensor as 



7-Oi 



J 
W 



W 



(24) 



(i.e., in terms of the comoving-frame radiation energy density 
J, momentum density W", and stress /C*-', respectively). In 
terms of the monochromatic energy density, momentum den- 
sity, and stress defined in Equation (24), the monochromatic 
lab-frame stress-energy tensor can now be written as 

T""" = L^g L\ J + {l\ L\ + L\ L^^jW + L'^-jlCl 

(25) 

We can also write JQ^ = fc'^ J, and define the monochromatic 
symmetric rank-two "variable Eddington tensor" as 



^ J,,n^n3fdn 



(26) 



where the monochromatic number-flux four-vector in the 
comoving-frame is defined in terms of angular moments of 
the distribution function as 



For the moment equations we derive in the next section it 
is useful to define the completely symmetric third-order mo- 
ment 



1 



(20) 



The solid angle element is dO, — sin -d d-d d(p, and the an- 
gular integration is carried out over the unit sphere; i.e., 

. . . dfl = J^^ In ■ ■ ■ d'd dip. The monochromatic 
lab-frame number-flux four-vector is related to the corre- 
sponding comoving-frame four-vector by 



(21) 



Jv^ \P\ 

= L'',L'',;L'', / U^'^Pe^de 



(27) 



(28) 



where the monochromatic third-order comoving-frame mo- 
ments are given by 



U^'''P = - / p'^p'^pPfdn, 



(29) 
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with U^'''P = L^^ L^f^ L'^-Ui^'^P. Note in particular that 

^bfL0 ^ ^A0£> ^ ui^M ^ £ q-fiP^ (30) 

and that we can expand the monochromatic third-order lab- 
frame moments in terms of the comoving-frame moments as 







LP,WL\LP.^ 



+L\ W i^g i^g + LP, W L'i L'-g + ^1 ^"j ^"o 

i'l IC '^ L^^ + L", LP- L'i + LP, L'^.LP-^ C^'' ,(31) 

where £*^'^ = Z*^'"' J', and we have introduced the completely 
symmetric rank-three variable Eddington tensor 



frfk ^ Jn^'n^n'' fdn 



(32) 



We have defined the lab-frame number-flux four-vector, the 
lab-frame stress-energy tensor, and the third-order lab-frame 
moments (a rank-three tensor). They are expressed in terms of 
angular moments of the distribution function in the comoving- 
frame (cf. Equations (21), (25), and (31)). From the defini- 
tions in Equations (20) and (23), we clearly have eAf^ ~ 

T'^°. This expression can be generalized to a useful covariant 
expression relating the number-flux four- vector and the stress- 
energy tensor by noting that T'^° = — WqT'^" = —UyTP'^ 

(Cardall et al. 2012). Thus, we have 



e 



(33) 



which is a covariant expression (valid in any frame; although 
e is always the neutrino energy measured by a comoving ob- 
server). Similarly, we have 



V 



1 



(34) 



which is also a covariant expression relating the stress-energy 
tensor and the third-order tensor 

4. GENERAL RELATIVISTIC ANGULAR MOMENT EQUATIONS 

We derive general relativistic angular moment equations in 
this section. The main results are the number-conservative 
monochromatic lab-frame number equation (Equation (41)), 
which is based on Equation (9), and the four-momentum- 
conservative monochromatic lab-frame stress-energy equa- 
tion (Equation (47)), which is based on Equation (12). The 
steps are straightforward. 

4.1. Number Equation 

From the definition of the monochromatic number-flux 
four- vector in Equations (20) and (21), and the expression in- 
side the spacetime derivative in Equation (9), it is clear that we 
obtain the evolution equation for the lab-frame neutrino num- 
ber density by integrating Equation (9) over the unit sphere Vl 
and dividing by e. For the spacetime divergence we find 



1 



-g dxP 



(35) 



With spherical momentum space coordinates, the momen- 



tum space divergence in Equation (9) becomes 



1 d 



P 9nJ ^ ^ ■> 



e sin d du 
1 d 



(36) 



We only focus on the first term (containing the derivative with 
respect to energy e) on the right-hand side of Equation (36). 
The second and third terms on the right-hand side of Equa- 
tion (36) vanish upon integration over momentum space an- 
gles {d,Lp). In particular, integrating over angles and dividing 
by e, the momentum space divergence becomes 



£2 de V 



ppp'-fdn], 



(37) 



where we have used the fact that 

e = Pi p^p^^e p^ p^ (38) 

(cf. Equation (16)). The rightmost expression in Equation 



(38) is derived from the mass-shell constraint (phpP 



0; 



i.e., p° is considered a function of the independent momen- 
tum space coordinates {p*}) and Equation (5); cf. Equation 
(136) in CM03. Then, writing the expression inside the en- 
ergy derivative in Equation (37) in terms of comoving-frame 
angular moments of the distribution function, we have 



2 T^O 



(39) 



where we have eliminated the connection coefficients associ- 
ated with the orthonormal comoving basis with Equation (7), 
and introduced the (coordinate basis) four-velocity of a co- 

moming observer u^, = L'^^uc, ~ ~L^^ (cf. Section 2). The 
covariant derivative of a vector Ai, is (cf. Landau & Lifshitz 
1975) 



dxP- 



VP A 



Thus, the number-conservative monochromatic lab-frame 
number equation, expressed explicitly in terms of angular mo- 
ments of the particle distribution function in the comoving- 
frame (via Equations (21) and (25)), becomes 



1 



d 



-J c[f] dn. 



Indeed, an integration over energy shells results in 



1 



d 



y/—g dxP- 
(cf. Lindquist 1966). 



C[/] 



d^p 
~\P\ 



(41) 



(42) 
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4.2. Stress-Energy Equation 

Similarly, from the definition of the stress-energy tensor in 
Equation (23) and the expression inside the spacetime deriva- 
tive in Equation (12), it is clear that we obtain the evolu- 
tion equation for the monochromatic lab-frame neutrino four- 
momentum by integrating Equation (12) over the full solid an- 
gle VL and dividing by e. The spacetime divergence becomes 



1 



d 



V~9 

The momentum space divergence in Equation (12) is 
1 d 



(43) 



e sin du^ 



" r M „p f 

i>P dpi ^ ^ ^ 



e de\ 

1 a 



siiii9a?9V "P dpi' ^ 



du 



.L^.p^p'pPf 



(44) 



Again, we focus only on the term with the energy derivative 
on the right-hand side of Equation (44). After integrating over 
the unit sphere and dividing by e, the momentum space diver- 
gence reduces to 



1 9 (.p po 



p^'p'^ppfdn 



(45) 



where the relations in Equation (38) have been used. We 
seek to express the momentum space divergence in terms of 
comoving-frame angular moments of the distribution func- 
tion, and write 



e L';, J^^ / p' pffdn^ ^ L^. r",^ 



(46) 



where we have taken steps similar to those used to arrive at 
Equation (39). 

Then, combining all the terms, the monochromatic 
lab-frame stress-energy equation, expressed in terms of 
comoving-frame angular moments of the distribution function 
(via Equations (25) and (31)), becomes 



pp -T-Pf 



- / 



[/] d^. (47) 



Equation (47) corresponds to Equation (3.18) in Shibata et al. 
(2011), which was derived from the moment formalism of 
Thome (1981). Integrating Equation (47) over comoving- 
frame energy shells gives the familiar result 



1 



d 



^T'-) +r^^,T'"' = / p^C [/] ^ (48) 



-g dx" 



(cf. Lindquist 1966). 

By invoking the transformation from the orthonormal tetrad 
basis to the coordinate basis we can rewrite Equation (47) 



in an equivalent form in terms of the stress-energy tensor as- 
sociated with the orthonormal tetrad basis 



yP q-pv 
^ pu ' 




P^C[/] dn, (49) 



where the stress-energy tensor T''" is related to the 
comoving-frame angular moments through the Lorentz trans- 
formation (Appendix A), the covariant derivative of a four- 
vector Ap is 



pP A_ 



(50) 



and the connection coefficients associated with the orthonor- 
mal tetrad basis are 



pp__ 



py 



. eP e-"- 



(51) 



Going even one step further, by invoking the Lorentz trans- 
formation, we can rewrite Equation (49) in an equivalent form 
in terms of the stress-energy tensor associated with the or- 
thonormal comoving basis 



1 



d 



pu 



yP rr-pv 
pv ' 



-g dx" 

-^l-^(e^U^'P^,Up) ^ - I pPC[f] dn, (52) 



where the covariant derivative of a four- vector is 



Vsv4p 



(53) 



and the connection coefficients associated with the orthonor- 
mal comoving basis are given in terms of connection coeffi- 
cients associated with the orthonormal tetrad basis by 



yP 

pi) 



u - V , (54) 

or in terms of the connection coefficients associated with the 
coordinate basis in Equation (7). 

The covariant nature of the general relativistic moment 
equations is illustrated by Equations (47), (49), and (52). Al- 
though they are analytically equivalent. Equation (47) may 
be preferred over Equations (49) and (52) as the basis for 
developing a numerical method to evolve the neutrino radia- 
tion field in full general relativity employing the so-called 3h-1 
decomposition (Shibata et al. 201 1; Cardall et al. 2012). In 
Section 5 we use Equation (47) to obtain moment equa- 
tions valid for conformally flat spacetimes. Equations (47) 
and (49) are conservative equations for the lab-frame four- 
momentum. In the absence of gravity and neutrino-matter 
interactions. Equations (47) and (49) express exact conser- 
vation of neutrino energy and — if Cartesian coordinates are 
used — momentum (Appendix A). Equation (52) is expressed 
in terms of the comoving-frame stress-energy tensor Re- 
cently, Miilleretal. (2010) presented numerical methods 
for neutrino radiation transport based on equations that are 
closely related to Equation (52) (see Appendix B). Equa- 
tion (52) is nonconservative — even in the absence of gravity 
and neutrino-matter interactions — since the connection coef- 
ficients depend on time and space derivatives of the fluid 
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three-velocity (cf. Equation (54)). In Appendix B we use 
Equation (47) to obtain conservative general relativistic mo- 
ment equations valid for spherically symmetric spacetimes 
(i.e., similar to Muller et al. 2010, but in conservative form). 

4.3. The Relationship Between the Number Equation and the 
Stress-Energy Equation 

In this subsection we illustrate the relationship between the 
moment equations for the neutrino four-momentum density 
and the moment equation for the neutrino number density; 
Equations (47) and (41), respectively. From Equation (33), 
it is evident how the neutrino number equation is related to 

, with 



the stress-energy equation. 
Equation (47) results in 



Thus, contracting — e 




C [/] dn. (55) 



To arrive at Equation (55), we have used Equations (33) and 
(34) in the first two terms on the left-hand side, and the co- 
variant relation u^,p^^ = u^p^ = — e to obtain the colli- 
sion term on the right-hand side. The third term on the left- 
hand side of Equation (55) consists of the 'leftover' terms af- 
ter bringing — e^^ u^j inside the spacetime and energy deriva- 
tives, and vanishes exactly by virtue of Equation (34). Thus, 
as expected. Equation (55) — the contraction of — e^^ with 
Equation (47) — reduces to Equation (41); i.e., the number- 
conservative monochromatic number equation. The solution 
to Equation (47) can thus be used to construct the solution to 
Equation (41) via Equation (33). 

When carrying out the steps to obtain the number- 
conservative monochromatic neutrino number equation from 
the four-momentum-conservative monochromatic stress- 
energy equation, we observe that terms emanating from the 
spacetime divergences and the geometry sources cancel with 
terms emanating from the momentum space divergences. The 
remaining terms constitute the left-hand side of the number 
equation. Such detailed knowledge about how the conserva- 
tive number equation is 'built into' the four-momentum con- 
servative stress-energy equation may be useful when develop- 
ing numerical methods. In particular, when discretizing the 
stress-energy equation to develop a two-moment model for 
neutrino transport, care should be taken to mimic the cancel- 
lations that occur in the continuum limit in order to ensure 
that the resulting numerical solution is also consistent with 
a discrete version of the conservative number equation, and 
thereby facilitating lepton number conservation as a 'built-in' 
property of the numerical method. We will elaborate further 
on this issue with specific examples in Sections 5 and 7, and 
in Appendix A. 

We have derived conservative general relativistic moment 
equations for the number density (Equation (41)) and 
the four-momentum density (e.g.. Equation (47)). The 
equations are in conservative form in the sense that (modulo 
neutrino-matter interactions and the geometry source terms in 
the stress-energy equation) the time rates of change of 
and T''" are governed by space and momentum space diver- 
gences. When integrated over comoving-frame energy bins 
(with de as the integration measure), the equations reduce 
to familiar position space conservation laws (Equations (42) 
and (48), respectively). The equations govern the evolution 
of quantities defined with respect to the global coordinate 



basis, which are explicitly expressed in terms of comoving- 
frame angular moments, M'^ and T'^'^ via Equations (21) and 
(25). Such equations may be suitable for implementation of 
neutrino transport capabilities in numerical codes intended 
for simulations of core-collapse supernovae, where accept- 
able conservation of total energy and lepton number is critical. 
We have also illustrated the relationship between the number 
equation and the stress-energy equation. 

In terms of a hierarchy of moment equations, we have 
truncated the hierarchy at the level of the so-called first- 
order moment (the comoving-frame energy flux density T-V"; 
the comoving-frame energy density J is proportional to the 
zeroth-order moment), which results in a two-moment model 
for neutrino transport. When truncating at this particular 
level, the resulting moment equations involve higher-order 
moments (or ratios of higher-order moments to the zeroth 
order moment; i.e., the variable Eddington tensors) fc*-' and 

r^^, which must be determined in order to obtain a closed 
system of equations. Due to symmetry, six unique com- 
ponents of fc'^ and ten unique components of must be 
determined (sixteen in total) per energy group. The pro- 
cedure to determine these higher-order moments in a clo- 
sure prescription that relates the higher moments to the first 
two is referred to as the moment closure problem. (We 
do not address the moment closure problem in this paper) 
Common approaches to the moment closure problem include 
analytic closure (e.g., Minerbo 1978; Levermore 1984), 
entropy-based closure (e.g., Cernohorsky & Bludman 1994; 
Levermore 1996; Smit et al. 2000; Brunner & Holloway 
2001; Hauck& McClarren 2010), and closure based on the 
solution of an approximate (simplified) Boltzmann equation 
(e.g., Rampp & Janka 2002; Miiller et al. 2010). 

5. MOMENT EQUATIONS FOR CONFORMALLY FLAT SPACETIMES 

The moment equations derived in Section 4 are valid for 
general spacetime coordinate systems. In this section we 
present moment equations suitable for simulations assum- 
ing a simplified spacetime metric. In particular, we impose 
the conformal flatness condition on the spatial metric (cf. 
Wilson etal. 1996; Isenberg 2008). The equations derived 
in this section are sufficiently general to accommodate the 
pseudo-Newtonian moment equations (Section 6), the New- 
tonian gravity, 0{v) approximation of the radiation moment 
equations (Section 7), the special relativistic limit (Appendix 
A), and the fully general relativistic case for spherically sym- 
metric spacetimes (Appendix B). (See Shibata et al. 2011; 
Cardall et al. 20 1 2, for moment equations for radiation trans- 
port employing the full 3h-1 formulation of general relativity.) 
To this end, we adopt the 3h-1 form of the spacetime metric 



ft- 
1^3 



(56) 



where a, and ^ij are the lapse function, the shift vec- 
tor, and the spatial three-metric, respectively. The confor- 
mal flatness condition is imposed by performing a confor- 
mal transformation of the three-metric, 7^- = "0 7y-, where 
?/'(.t'^) is the conformal factor, and insisting that the con- 
formally related metric 7^^ is diagonal. In particular, we 
set 7ij = diag [l, a^(x^), c^(a;^) ], where the metric 

functions (or scale factors) a, h, and c are sufficiently general 
to accommodate Cartesian, spherical, or cylindrical coordi- 
nates (see for example Table I in Cardall et al. 2005). The 
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scale factors, and hence the conformally related metric, may 
depend on the spatial position coordinate, but are taken here 
to be independent of time. The lapse function, the shift vector, 
and the conformal factor completely determine the spacetime 
metric, and can be obtained by solving a system of nonlinear 
elliptic equations (Wilson et al. 1996). The four-metric and 
its inverse lower and raise indices on four-vectors and tensors, 
while the spatial metric and its inverse 7'-' can be used to 
lower and raise spatial indices of purely spatial vectors and 
tensors. The determinant of the spacetime metric is denoted 
g = —o? 7, where 7 = V'^^ (a6c)^ is the determinant of the 
spatial metric. 

In the 3h-1 formulation of general relativity, the four- 
dimensional spacetime is sliced into a "stack" of spatial 
(spacelike) hypersurfaces of constant (global) time coordinate 
(see for example the recent book by Baumgarte & Shapiro 
2010). In particular, the unit normal to the spacelike hyper- 
surfaces is denoted = ( , ~a~^ ). With the metric 
specified in Equation (56) we have = g^u n'' = ( —a, ). 
Observers whose four-velocity coincides with 71^ are at rest 
with respect to the spacelike slice, and are referred to as Eule- 
rian observers (e.g., Banyuls et al. 1997). For the transforma- 
tion from the orthonormal tetrad basis to the global coordinate 
basis we can use 



a 







(57) 



where e\ = V ^diag [l, a ^, ^] (i.e., e'^e'^p.g^,^ = 

'7pp)- With e'^ given by Equation (57), we have e'^ = — 1. 
That is, the spacetime orientation of the orthonormal tetrad 
basis has been chosen so that e'g is aligned with the unit nor- 
mal of the spacelike hypersurfaces (i.e., ~ n^^). The four- 
velocity of a comoving observer with respect to the lab-frame 
coodinate basis can now be computed. The result is 



(58) 



where w'' = (O, w'j and = e\w' is the coordinate ba- 
sis three-velocity of the observer comoving with the fluid 
(cf. Appendix A). Note that = —W"^ ( 1 — ) 

^ (1 - v,v'y\ 

Foflowing Cardalletal. (2012), we form 'Eulerian' and 
'Lagrangian' decompositions of four-vectors and tensors. In 
particular, the Eulerian and Lagrangian decompositions of the 
coordinate basis number-flux four-vector are 



(59) 
(60) 



where the lab-frame (or Eulerian-frame) number density and 
flux (the 'Eulerian projections' of the number-flux four- 
vector) are — and Fj^ = e'^A/^, respectively. The 
comoving-frame number density and flux (the 'Lagrangian 

projections' of the number-flux four-vector) are Jj^ = 
and -Hj!^ = i^j7V'~ (cf. Equation (21)). Note that J"/ is or- 
thogonal to the Eulerian observer's four- velocity {n^Fj^ = 
0), while Hjj^ is orthogonal to the four- velocity of the comov- 
ing observer (u^ 'Hjv = 0). 

Similarly, the Eulerian and Lagrangian decompositions of 
the coordinate basis stress-energy tensor are 



(61) 
(62) 



where the monochromatic lab-frame radiation energy density, 
flux, and stress (the Eulerian projections of the stress-energy 
tensor) are related to the components of the stress-energy ten- 
sor in the orthonormal tetrad basis by 



T^^T" = e'jT*", and S^"" 



(63) 



respectively. Note that ri^ F'^ = S^'^ = S^'^ = 0. 
The monochromatic comoving-frame radiation energy den- 
sity, flux, and stress (the Lagrangian projections of the stress- 
energy tensor) are related to the components of the stress- 
energy tensor in the orthonormal comoving basis by (cf. 
Equation (25)) 



J = r™ , = L^j and /C^" = L'l L"^ V\ (64) 



where W = /C^'^ = IC" = 0. 

Finally, the Eulerian and Lagrangian decompositions of the 
third-order moment (rank-three tensor) are 



Ut'^P = g n'' n" nP + nP + V nP + TP n^" v}' 

nP + VP n"" + VP + Qp^p (65) 

uP + JQt'P u" + JCP uP + Ct'^P ) , (66) 



where g = U°°°, ^ e^\W°°, VP"" = e'^e^jW^, 
QP-'P = e'\e''^e\U'^'^\ a.nA Cp^p = Lf'-L^-LP-^ C^'' . Again, 
I'', VP", and Q^'^'^P are purely spatial and orthogonal to the 
four-velocity of the Eulerian observer, while CP'^p is orthog- 
onal to the four-velocity of the comoving observer The 
Eulerian and Lagrangian decompositions of the number-flux 
four-vector, the radiation stress-energy tensor, and the third- 
order moments are useful when deriving conservative moment 
equations for the radiation field — in particular, when seek- 
ing conservative evolution equations for lab-frame radiation 
quantities, expressed in terms of comoving-frame moments. 

The radiation moment equations for conformally flat space- 
times can be obtained by specializing the 3h-1 moment equa- 
tions in Cardall et al. (2012), or by adopting the spacetime 
metric in Equation (56), with jij = ^j'^jij, and computing 
connection coefficients directly. Here we have adopted the 
latter approach, but see also Cardall et al. (2012). We obtain 
the evolution equation for the monochromatic lab-frame radi- 
ation energy density by contracting Equation (47) with — n^. 
The result is 



1 d 



'gdt 



-g dx^ 



9 In a 



S' 



d In ijp' 
' dr 



1 <9 / 9 



(67) 
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where we have defined the energy space energy flux 



dv'- 



dr 

V - Qv 



-pn 
a 



9 In a 

9 In f/;^ 



dW 
d'r 



ft 



'dlnih'^ ,; d In ^|? 

1— + yJ _ 

^ dr dx^ 

+ in,. 



a dx 
dW 

9^ 



-I 



(68) 



We have also defined the covariant derivative with respect to 
the conformally related three-metric 



where 



1 

9j 



kl 



dx^ 



(69) 



(70) 



dx"^ dx'- 

and the proper time derivative along constant coordinate lines 

d Id 13' d 



dr 



dt 



da 



(71) 



We obtain the monochromatic radiation momentum equa- 
tion by contracting Equation (47) with the projection operator 
7v = 5v + "m- The result is 



1 d 
7=ldi 



1 



d 



^ 2 dx^ 
^_d_ 

de 



£ 



sj —g dx^ 
9 In a 



x/7 



dx' 



J d In ^jj"^ 



1 dP^ 



1 



dx' 

pp c [/] dn, 



dx' 



(72) 



where we have defined the energy space momentum flux 

r„ <9u^ ^ I, dv^ 1 / ^ 



dx^ 
9 In a 
dx^ 
^ 9 In '0'^ 
dx^ 
\dW 

V a7 



' (9 In i\p- 



, 9 In 



9t 



(73) 



Equations (67) and (72) are conservative equations for the lab- 
frame radiation energy density and momentum density. The 
independent variables are the coordinate basis spacetime po- 
sition components x^ and the neutrino energy measured by 
a comoving observer e. When integrated over the comoving- 
frame energy, the equations reduce to familiar position space 
conservation laws. They express exact energy and momentum 
conservation in the absence of neutrino-matter interactions 
and geometry sources (due to gravity and curvilinear coordi- 
nates). The terms in the energy derivatives (cf. Equations (68) 
and (73)) result in changes in the radiation energy spectrum 
measured in the comoving-frame from gravitational redshifts 



and acceleration of the comoving observer Note that Equa- 
tions (67) and (72) differ from the 3h-1 general relativistic mo- 
ment equations given by Shibata et al. (201 1); Cardall et al. 
(2012), in that they are not expressed in terms of the extrinsic 
curvature Kij (an evolved quantity in 3h-1 general relativity; 
e.g., Baumgarte & Shapiro 2010). 

For the terms inside the energy derivatives in Equations 
(67) and (72), we can use Equation (34), and write = 
ly ( 71^ + ), to relate Eulerian projections of the third- 
order moments (Equation (65)) in terms of Eulerian projec- 
tions of the stress-energy tensor (Equation (61); Cardall et al. 
2012). In particular, we have 



WV - e| + S'^v, + — Q'J^,ji;fe|, (74) 

WV' = t\s'' + ■^Q'^'^i^fc}, (75) 
W{V - gu') = e{(j-' - fw') + [S'l - vT'^Vj 
^{^Q^o^ „ ^^S^A^^.Vk - v'^Q^'^'v.vuviY (76) 

+ 2*'" - - Q'^'vkvi } , (77) 

w(i',T -G) ^ (78) 

wi^VjV'^ -r^ = -er, (79) 

W (ufc Q'^'' - P*^ ) = - e 5'^ . (80) 



The neutrino number equation can be obtained from Equa- 
tions (67) and (72). From Equation (33) we have 



W / 



(81) 
(82) 



As was done in Section 4.3, we obtain a conservative equation 
for the monochromatic lab-frame neutrino number density by 
adding (T^W times Equation (67) and —eWv' contracted 
with Equation (72) — and subsequently bringing the necessary 
terms inside the time, space, and energy derivatives, respec- 
tively. The result is 



1 



V7 



V-gdx 
-I C [/] 



(83) 
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where the energy space number flux is 

W 



-W 



9 In a /dlnip'^ 



1 



.9% 



5^ 



9 In tp^ 
dW 



S' 



dx^ 



1 9/3- 



P - v'S. 



^ a dx 
) dx3 



-I 



(84) 



When deriving Equation (83) from Equations (67) and (72), 
terms emanating from the time derivatives, space derivatives, 
and the geometry sources cancel exactly with terms emanat- 
ing from the energy derivatives (see also Cardall et al. 2012). 
Ideally, a discrete representation of the neutrino energy and 
momentum equations can be constructed so that a conserva- 
tive neutrino number equation can be analogously obtained in 
the discrete limit. The discretization of the left-hand sides of 
the energy and momentum equations is then consistent with 
neutrino number conservation. (A similar consistency must 
also be considered for the right-hand side ot the equations.) 
Clearly, in order to construct such a discretization, the indi- 
vidual terms in the neutrino energy and momentum equations 
cannot be discretized independently. 

6. PSEUDO-NEWTONIAN MOMENT EQUATIONS FOR NEUTRINO 
RADIATION TRANSPORT 

In this section, we specialize the moment equations pre- 
sented in Section 5 to the 'pseudo-Newtonian' limit by adopt- 
ing a spacetime metric consistent with the line element (cf. 
Kimetal. 2009,2012) 

ds^ = - ( 1 + 2 $ ) <i<2 + ( 1 + 2 $ dx' dx\ (85) 

where $ <C 1 is a pseudo-Newtonian gravitational po- 
tential. In this approximation, all orders in the fluid ve- 
locity V are retained, and the pseudo-Newtonian gravita- 
tional potential is obtained by solving a Poisson equation 
in which the trace of the stress-energy tensor — not just the 
rest-mass density — contributes to the source on the right- 
hand side. We allow for common curvilinear spatial coor- 
dinates (i.e., Cartesian, spherical, and cylidrical) by setting 
% = diag[l, a^{x'^), b"^ {x'^) {x^)] (Section 5). The spa- 
tial metric 7^ = (1 + 2 4>)~^ and its inverse are used to 
lower and raise indices on (spatial) vectors and tensors. The 
determinant of the spacetime metric becomes 



(l + 2<I>)-iy^. 



(86) 



Similarly, we have = (1 + 2 <I>)^'^/^-y/7, where = 
abc. Note that the spacetime metric in Equation (85) dif- 
fers slightly from the Newtonian one given by Misner et al. 
(1973);Schutz (1985), but agrees to £)($); i.e., (l+2$)-i = 
(1 — 2$) + 0($^). Using the spacetime metric in Equation 
(85) turns out to be algebraically advantageous. Note in par- 
ticular that ^/~^ = <5 V7> where a = (1 + 2$)^/^. We 
note in passing that Kim et al. (2009, 2012) recently com- 
bined relativistic hydrodynamics with the weak-field (Newto- 
nian) limit of general relativity (cf. Equation (85)) to study 
equilibrium solutions of rotating relativistic stars and found 
good agreement with general relativistic computations in the 



reported cases. (See also Takiwaki et al. 2009, who adopted 
a similar approach in magneto-rotational core-collapse simu- 
lations.) 

We obtain the pseudo-Newtonian radiation moment equa- 
tions from Equations (67) and (72) by setting a = a = 
(1 + 2<I>)i/2, = 0, and i})^ = = (1 + 2<^)-^/'^, and 
retaining all terms of 0{^). We state the results for easy ref- 
erence. The radiation energy equation becomes 



1 d 



~9dt 



1 



d 



-g dx^ 



(9$ 



-gr — 



dx^ 



where the lab-frame radiation energy density, momentum 
density, and stress {£, T', and iS'-* , respectively) are related to 
the comoving-frame moments in Equations (A9)-(A11). The 
energy space energy flux (cf. Equation (68)) reduces to 



. 9% 
dx'' 



dx^. 
dW 



Similarly, the radiation momentum equation becomes 



(88) 



1 d 
7=gdt 



1 



d 



dx^ 



2 



,9$ 

dx^ 



^ dx^ dt 



7^,- / yC[/] dO, 
e Jn 



(89) 



where the energy space momentum flux is given by (cf. Equa- 
tion (73)) 



dx'- 

O . I 1- ir — - 



dW 



Equations (87) and (89) are vaUd to all orders in v, but hm- 
ited to weak gravitational fields ($ ^ 1). In the case of no 
gravitational fields ($ = 0) they reduce to the special rela- 
tivistic moment equations (cf. Equations (A22) and (A24) in 
Appendix A). Moreover, Equations (87) and (89) are consis- 
tent with the conservative neutrino number equation given by 



1 



1 



d 



V-gdx 



(91) 
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where the energy space number flux is (cf. Equation (84)) 



The lab-frame number density and number flux are related 
to the corresponding comoving-frame moments in Equations 
(A6) and (A7), respectively. 

Equations (87), (89), and (91) simplify even further when 
slowly varying gravitational fields are considered (d^/dt = 
0;cf. Kimetal. 2012). 

7. NON-RELATIVISTIC SELF-GRAVITATING NEUTRINO 
RADIATION HYDRODYNAMICS 

In this section we detail the full set of neutrino radia- 
tion hydrodynamics equations intended for our planned non- 
relativistic simulations of core-collapse supemovae and re- 
lated systems. The equations are deduced from the relativistic 
equations derived in previous sections. Per our initial dis- 
cussion in Section 1, in order to formulate a system of mo- 
ment equations that are consistent with a conservative equa- 
tion for the neutrino number density, we adopt different orders 
of V for the radiation energy and momentum equations. We 
consider two possible cases: the 'C'(v)-plus' and the '0{v)- 
minus' moment equations. The 0{v)-plus moment equations 
consist of an 0{v'^) energy equation and an 0{v) momentum 
equation, and are consistent with the conservative O(u^) neu- 
trino number equation. Similarly, the C'(w)-minus moment 
equations consist of an 0{v) energy equation and an 0(1) 
momentum equation, and are consistent with the conservative 
0{v) neutrino number equation. Dimensional analysis sug- 
gests that gravitational effects on the radiation field are for- 
mally 0{v'^). Therefore, in the C'(v)-plus moment equations 
for neutrino radiation transport, we also retain terms due to 
a Newtonian gravitational potential. These 'gravitational red- 
shift' terms may play a non-negligible role in the post-bounce 
supernova environment (e.g., Bruenn et al. 2001). More- 
over, by including gravitational redshift effects, the 0{v)- 
plus moment equations are conceptually more similar to the 
fully relativistic moment equations (e.g., Shibata et al. 201 1; 
Cardalletal. 2012). 

7.1. 0{v)-Plus Moment Equations for Neutrino Transport 

In this subsection we detail radiation moment equations for- 
mally valid to 0{v), including effects due to a Newtonian 
gravitational field. We obtain these Newtonian, C'(u)-plus ra- 
diation moment equations from the pseudo-Newtonian equa- 
tions in Section 6 — Equations (87) and (89) — ^by considering 
non-relativistic fluid velocities. In the radiation moment equa- 
tions presented in this section, we also omit terms contain- 
ing the time derivative of the gravitational potential. Further- 
more, we retain at least all terms that are linear in the fluid 
velocity. We retain O(w^) terms in the radiation energy equa- 
tion, while we retain 0{v) terms in the radiation momentum 
equation. As discussed in the introduction, and in detail be- 
low, we adopt this 'mixed' order of the radiation energy and 
momentum equations in order to achieve satisfactory consis- 
tency with the conservative neutrino number equation. Since 



we assume v <^ 1, we set W 



W 



1 + , where 



= v''Vi. Then, the monochromatic O(u^) radiation energy 
equation becomes (cf. Equation (87)) 



(93) 



where the 0{v'^) lab-frame radiation energy density and mo- 
mentum density are related to the comoving-frame moments 
{J, W, and JC'^) by 



(94) 
(95) 



and the 0{v) radiation energy density, momentum density, 
and stress are related to the comoving-frame moments by 



£ = J + 2v,n\ 



(96) 
(97) 
(98) 



(cf. Equations (A9)-(A11)). Radiation quantities adorned 
with a tilde (e.g., £) are accurate to 0{v'^). To 0{v) we also 
have 



(99) 
(100) 

(101) 



(cf. Equations (A12), (A13), and (A15)). The comoving- 
frame moments J, W, JC^^, and are defined in Section 3. 
The energy space energy flux, written in terms of comoving- 
frame moments, is (cf. Equation (88)) 



2v'IC,j]- — 



1 dv^ 
a dt 



V, + Hj v'' + v' C " 



k , ..i rjk\jd2jk 
■ \ 0$ 



~ld\nW ~:dhiW 

-£ j:^ 

a dt dx'- 



}■ 



(102) 



where the 0{v) and 0{^) terms are defined separately in 

dx3 /■ 



(103) 

In Equation (102), except for the last two terms, we have 
dropped terms of order equal to or higher than 0{v^) and 
0{v'^ $). Similarly, the monochromatic 0{v) radiation mo- 
mentum equation becomes 



1 ^ ( TT 



2 dx^ 
C^, (104) 
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where energy space momentum flux is 



L a 



1 9 In 



dlnW 



dt ' dxi 

and the 0{v) and 0($) terms are contained in 



. W -I 



(105) 



L a 



dt dx^^2^' " dx^ " 



dxir 

(106) 

In Equation (105), except for the last two terms, we have 
dropped terms of order equal to or higher than 0{v'^) and 
0{v Note that 5'^ is not accurate to 0{v^), but is adorned 
with a tilde to signify that some higher-order terms have been 
retained. Also note that we have retained the ^-dependence 
in the determinant of the spacetime metric (cf. Equation (86)) 
appearing in the moment equations, and retained factors a and 
tjj, which give rise to higher-order terms in some cases (e.g.. 
Equation (106)). These should all be kept for consistency with 
the conservative 0{v'^) number equation (see below). 

In Equations (93) and (104), the collision terms on the right- 
hand side are, to 0{v'^) and 0{v), respectively. 



In Jn 



[/] dn, 



(107) 



2^C[f]dn. (108) 



(vi is the fluid three-velocity with respect to the orthonormal 
tetrad basis; Appendix A.) 

Equations (93) and (104) are conservative equations for the 
monochromatic lab-frame radiation energy density and mo- 
mentum density, respectively. Equation (93) expresses ex- 
act conservation of radiation energy in the absence of gravity 
and neutrino-matter interactions, and — if in addition, Carte- 
sian coordinates are used — Equation (104) expresses exact 
conservation of radiation momentum. In the strict 0{v) 
limit'*. Equations (93) and (104) can be compared with corre- 
sponding equations derived by Buchler (1979, his Equations 
(9) and (10), respectively). See also Kanekoetal. (1984); 
Munier & Weaver (1986b); Buras et al. (2006). We obtain 
the strict 0{v) limit of Equations (93) and (104) by set- 
ting $ = 0, letting i,P,P,S\ £,J'\J"',S% and 
C'^ by setting W = I in Equation (107). Buchler's 

equations evolve the comoving-frame moments and are non- 
conservative. However, to 0{v), the terms inside the energy 
derivatives in Equations (93) and (104) agree with the energy 
derivative terms in Buchler's energy and momentum equa- 
tions. Shibataetal. (2011) also list conservative moment 
equations for neutrino radiation transport in the 'slow-motion' 
limit (cf. their Equations (8.5) and (8.6)). They do not include 
gravitational effects, or any higher-order terms in the fluid 
velocity. They also do not include the terms involving the 
derivative of the fluid velocity with respect to time in the en- 
ergy derivative terms, which can be important when A/r < 1 
(e.g., in optically thin regions; Buchler 1979), where A and r 
are typical length and time scales, respectively. 

By integrating Equations (93) and (104) over the comoving- 
frame energy (with de as the measure of integration) we 

By the "strict" 0{v) limit we mean that all (or most) tenns that are at 
most linear in the fluid velocity have been retained in the radiation energy 
and momentum equations. As we emphasize in this section, in this strict 
0(v) limit, the radiation moment equations are not fully consistent with a 
consen'citive equation for the neutrino number density. 



obtain the energy-integrated (grey) Eulerian-frame radiation 
energy equation. 



9$ 

ox^ 




(109) 



and the grey Eulerian-frame radiation momentum equation. 



1 di^FA 1 di^gS- 



dt 



-g dx-> 



\ 9$ 



where the grey radiation variables are 



2 dx^ 
(110) 



{E,F\F\ S'^ , J, K'^ , , } 

J-' , S'^ , J, , ,C^'}e^de. (Ill) 



Note that we use the font type to distinguish the energy- 
dependent radiation variables (denoted with calligraphic font) 
from the grey radiation variables. The Eulerian-frame grey 
radiation energy density, momentum density, and stress {E, 
F\ and S"*^, respectively) are related to the corresponding 
comoving-frame quantities (J, H^, and A'*-') via relations 
analogous to those listed in Equations (96)-(98). In the strict 
0{v) limit, for Cartesian coordinates, and with $ = 0, Equa- 
tions (109) and (110) correspond to Equations (32a) and (32b) 
given by Lowrie et al. (2001), which they refer to as the 'cor- 
rect' 0{v) Eulerian-frame radiation energy and momentum 
equations, respectively. One requirement for the correct 0{v) 
limit is that the hyperbolic wave speeds associated with the 
system of equations is bounded by the speed of light. The 
presence of geometry sources due to the gravitational field 
or the use of curvilinear coordinates does not destroy the hy- 
perbolic character of the equations, since they do not con- 
tain any differential operators acting on the radiation variables 
(Banyuls et al. 1997). 

Equations (93) and (104) are correct to 0{v'^) and 0{v), 
respectively. We have explicitly retained higher-order (non- 
linear) terms in the fluid velocity inside the energy deriva- 
tive of both equations. In the energy equation, we have also 
retained terms containing the dot product of the fluid three- 
velocity with the gradient of the gravitational potential. How- 
ever, when combined, the system is only accurate to 0{v). 
Hence, we refer to these equations as the '0(i))-plus' approx- 
imation of the radiation moment equations. These higher- 
order terms must (by definition) be small in the 0{v) limit 
we are interested in here. However, by retaining these terms, 
the solution to Equations (93) and (104) becomes consistent 
with the conservative O(w^) neutrino number equation. This 
consistency may help enable exact lepton number conserva- 
tion in the two-moment model for neutrino radiation transport 
in the 0{v) limit. We elaborate further on the details here. 

Equation (81) relates the Eulerian-frame neutrino number 
density to the corresponding neutrino energy density and mo- 
mentum density. The equations governing their time evolu- 
tion are similarly related. In particular, by adding e^^W times 
Equation (93) and —e~^Wv'^ contracted with Equation (104) 
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we obtain 



1 d 




W_ 

e 



a ot 



from the time derivative terms. 
Similarly we obtain 



J a at 
Remember W = 



(112) 




dW 



(113) 



from the space derivative terms, while the combination of the 
geometry sources results in 



W f 1 r 



+ 



(h^+z;'/C/+/C>^)0}. (114) 



From the energy derivative terms we obtain 
1_ d 

where the number-flux in energy space is obtained from 



-2 5e(' 



2 P- 

■J K 



(115) 



W 



w 



f Jv, - 



1 dv^ 
a dt 



J dxi 



dx^ 



/ dx^ / 

Finally, the collision terms combine to give 



dW 



W 



= i- 



(117) 



(Note that Vi e\ (5\ = vf. Appendix A.) 

When adding Equations (112)-(115), which result in the 
left-hand side of the conservative neutrino number equation 
(Equation (1 18) below), we note that the second term in Equa- 
tion (115) (cf. Equation (116)) cancels exactly with the 'left- 
over' terms in Equations (112) and (113), and Equation (1 14). 
In particular, the first term in the second line of Equation (112) 
cancels with terms emanating from the second line on the 
right-hand side of Equation (116). Similarly, the first term 
in the second line of Equation (113) cancels with terms em- 
anating from the third line on the right-hand side of Equa- 
tion (116). The terms emanating from the geometry sources 
(Equation (1 14)) cancel with terms emanating from the fourth 
and fifth line on the right-hand side of Equation (116). Finally, 
the terms involving time and space derivatives of the approxi- 
mate Lorentz factor W cancel with the terms emanating from 



the last line on the right-hand side of Equation (116). More- 
over, when contracted with —We~^ u\ several terms in the 
energy derivative of the radiation momentum equation cancel 
exactly with some of the 0{v'^) terms in the energy derivative 
of the iye~^-multiplied radiation energy equation (these can- 
cellations are only obtained by retaining higher-order terms 
in Equation (93)). We also find that the 0{v'^) terms involv- 
ing the time and space derivatives of W in Equations (102) 
and (105) are needed to cancel with the corresponding terms 
in Equations (112) and (113), obtained after pulling W inside 
the time and space derivatives. 

By combining all the terms (Equations (1 12)-(1 17)) we ob- 
tain 



1 d 



d 



g dt\ I \/—g dx"^ 

£2 de \ •^^)^^ ' 



(118) 



where the Eulerian-frame 0{v'^) number density and number 
flux density are (cf. Equations (A6) and (A7)) 

£j^^—[£^ = WJj^ + v,h;^ + 0(w3),(119) 

= ^{^'- =U^ + v'Jm + O(f'), (120) 

where we have used Equations (94)-(95) and (96)-(98). We 
have nil — e\ (5\ "H and the comoving-frame moments J^^ 
and n I are defined in Section 3. We have retained the 0{v^) 
term appearing in the collision term on the right-hand side of 
Equation (118) (defined in Equation (117)), but this term can 
safely be dropped in practical computations. 

Equation (118) is a conservative equation for the Eulerian- 
frame neutrino number density. It is valid to 0{v'^ ). Together 
with an equation for the electron density, it states that the total 
lepton number is conserved during lepton number exchange 
with the fluid (see discussion in Section 7.4). In the absence 
of neutrino-matter interactions. Equation (118) states conser- 
vation of particle number. It is obtained analytically from the 
monochromatic radiation energy and momentum equations. 
Modulo terms of 0{v^) or higher (cf. Equations (1 17), (1 19), 
and (120)), it can also be obtained directly from Equation 
(83). 

A numerical solution to Equations (93) and (104) should 
also be consistent with Equation (118) in order to ensure 
lepton number conservation in simulations of neutrino radi- 
ation transport in core-collapse supernovae and related sys- 
tems (Section 7.4). Ideally, the discretized neutrino energy 
and momentum equations can be similarly combined to obtain 
a discrete representation of the conservative neutrino number 
equation. Note that we arrived at Equation (118) in an an- 
alytically exact manner (due to exact cancellations; we did 
not throw away any 0{v^) terms outside the time, space, or 
energy derivatives). By retaining the higher-order terms in 
Equation (93) and Equation (104), the two-moment model 
of neutrino radiation transport is consistent with the conser- 
vative 0(^2 ) neutrino number equation. If the higher-order 
terms in Equations (93) and (104) are not retained (i.e., in the 
strict 0{v) limit), there will be additional 0{v'^) terms in the 
number equation derived from the 0{v) energy and momen- 
tum equations, and the consistency of the two-moment model 
with the number equation is only 0{v) — which may be ac- 
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ceptable in practical numerical computations. However, note 
that by omitting the dv'^ / dt-teim inside the energy derivative 
in the energy equation, the consistency with the conservative 
neutrino number equation can reduce to 0(1). Thus, omit- 
ting this term in numerical simulations of core-collapse super- 
novae can potentially result in severe consequences for lepton 
number conservation. 

By integrating Equation (118) over the comoving-frame 
energy we obtain the Eulerian-frame grey neutrino number 
equation 



.gF; (121) 



where 



{ F^, } 



}e'^de. (122) 



We have presented conservative equations for the 
monochromatic lab-frame radiation energy density and mo- 
mentum density — valid to 0{v'^) and 0{v), respectively — for 
common curvilinear coordinates (Cartesian, spherical polar, 
and cylindrical); Equations (93) and (104), respectively. We 
have also demonstrated how these equations (and therefore 
also their solutions) are fully consistent with the conservative 
equation for the monochromatic 0{v'^) lab-frame radiation 
particle density (Equation (118)). Note that for a consistent 
description of neutrino radiation hydrodynamics with the 
C'(w)-plus moment equations, the hydrodynamics equations 
may also have to be promoted to include higher-order terms 
in the fluid velocity. 



7.2. Further Simplifications of the Neutrino Radiation 
Transport Equations: 0{v)-Minus Moment Equations 

In this subsection we further specialize the radiation mo- 
ment equations presented in Section 7.1 by introducing 
the C'(f)-minus moment equations. Although apparently 
less complex than the fully relativistic moment equations 
(cf. Shibataetal. 2011; Cardall et al. 2012, Section 5), 
the C'(t;)-plus moment equations presented in Section 7.1 
are still nontrivial to discretize for numerical computations. 
(Their complexity rivals that of the pseudo-Newtonian mo- 
ment equations presented in Section 6.) Therefore, we pro- 
pose a further simplification, which involves solving the 0{v) 
energy equation and the 0(1) momentum equation, as a use- 
ful intermediate (first) step beyond the 0(1) moment for- 
malism previously used by some to model neutrino radiation 
transport in core-collapse supernovae (e.g.. Burrows et al. 
2006). The resulting system is formally valid to 0(1). How- 
ever, the moment equations are consistent with the conser- 
vative 0{v) neutrino number equation. Moreover, the equa- 
tions presented here evolve the radiation energy density and 
momentum density, and, with a proper closure prescrip- 
tion, may be an improvement over the equations solved in 
common 0{v) multigroup flux-limited diffusion approaches 
(e.g., Swesty & Myra 2009; Zhang et al. 2012). In the non- 
relativistic, Newtonian gravity limit, we have <&, c^, ^ 1, 
where Cs is the sound speed. Thus, we ignore gravitational 
effects on the radiation field in this subsection (i.e., $ = 0). 

Then, the monochromatic 0{v) lab-frame radiation energy 



equation becomes (cf. Equation (93)) 



:[/]d^^ = c^ (123) 



where the 0{v) lab-frame radiation energy density and mo- 
mentum density, £ and T\ are given by Equations (96) and 
(97), respectively. The energy space energy flux F^ is given 
by Equation (103), with $ = 0. 

Similarly, the monochromatic 0(1) radiation momentum 
equation becomes (cf. Equation (104)) 



2 dx 



'jk 



(124) 



In Equation (124), the energy derivative terms vanish since 
we have dropped all the velocity-dependent terms in addition 
to the gravitational terms (cf. Equation (105)). 

Equations (123) and (124) are consistent with the conser- 
vative 0{v) neutrino number equation: by adding e^^ times 
Equation (123) and — e^^v' contracted with Equation (124) 
we obtain 

d£^r 1 d f ^ A Id 



dt yj^ dx 



de 



(125) 

where the 0{v) lab-frame number density and number flux 
density, and the number flux in energy space are 



2 dx^ 



„fc 



(126) 
(127) 

(128) 



Equation (125) is the conservative 0{v) number equation, 
and is obtained in a manner similar to the O (u^ ) number equa- 
tion detailed in Section 7.1. On the left-hand side, the extra 
terms emanating from pulling the fluid three-velocity inside 
the time and space derivatives of the momentum equation can- 
cel with the extra terms emanating from pulling e^^ inside the 
energy derivative of the energy equation. On the right-hand 
side, the velocity-dependent term in the collision term of the 
energy equation cancels with the contraction of — e^^u* with 
the collision term of the momentum equation. However, by 
reducing the order of the energy and momentum equations 
(to 0{v) and 0(1), respectively) the number of cancellations 
that occur is dramatically reduced. Because of these simplifi- 
cations. Equations (123) and (124), which are consistent with 
the conservative 0{v) neutrino number equation, may be a 
suitable starting point for developing lepton number conserva- 
tive numerical methods for neutrino radiation hydrodynamics 
based on the two-moment model. 

7.3. Non-Relativistic, Self-Gravitating Hydrodynamics 

For self-consistent Newtonian, non-relativistic simulations 
of self-gravitating neutrino radiation hydrodynamics, the mo- 
ment equations for the radiation field in Section 7.2 (one set 
for each of the neutrino species) must be coupled to the equa- 
tions of self-gravitating hydrodynamics. For completeness 
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we list the non-relativistic hydrodynamics equations includ- 
ing self-gravity in this subsection. (For brevity do we not 
consider nuclear reactions here. However, see for example 
Plewa & Miiller 1999.) We consider a perfect fluid; i.e., we 
ignore fluid viscosity and thermal conduction. 

The equations describing a self-gravitating perfect fluid 
(e.g.. Landau & Lifshitz 1959) include the mass conservation 
equation 



dp I d , , 



0, 



the fluid momentum equation 



dt dxi 
and the fluid energy equation 



2 dx' 



(129) 



(130) 



J__d_ 

dt -^7 dx^ 



-V 



-pv 



dx^ 



C. (131) 



For simulations involving a nuclear equation of state. Equa- 
tions (129)-(131) must be supplied with a balance equation 
for the electron number density 



drie 1 9 / _ 



TieV = 



= - C" - CP. 



(132) 



In Equations (129)-(132), the fluid energy density, momen- 
tum density and stress are E^ = e + ^pviv'^, F^i = pvi, 
and j — pv^ Vj + p, where p — frib rib, "w', e, and rig 
denote the mass density, the ith component of the fluid three- 
velocity, the fluid pressure and internal energy density, and 
the electron density (electrons minus positrons), respectively. 
The average baryon mass and the baryon density are denoted 
TTib and rife. Equations (129)-(132) are closed with the spec- 
ification of an equation of state p ~ p{p,T,Ye), where T 
and Ye = ne/ni, are the fluid temperature and electron frac- 
tion, respectively. In Equations (130) and (131), the collision 
terms include energy and momentum exchange with all neu- 
trino species s; i.e.. 



{c-,c-} = 5]{cf,c-}. 



(133) 



Here we consider electron lepton number exchange between 
the fluid and the neutrino radiation field. Only interactions 
involving electron neutrinos (i^e) and electron antineutrinos 
(De) contribute to the right-hand side of Equation (132). 

The Newtonian gravitational potential is obtained by solv- 
ing Poisson's equation 



^/j dx^ \ dx^ / 



Att p. 



(134) 



Note that in the non-relativistic, Newtonian gravity limit only 
the mass density contributes to the gravitational field (as op- 
posed to all types of stress-energy in general relativity); i.e., 
we assume w <C 1 and e,p, J, K\ ^ p. 

Using Equation (134), we can rewrite the gravitational force 
in Equation (130) as 
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1 q ,fc dljk 
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(135) 



where the gravitational stress tensor is defined as 



^4, J - 



47r V 



(136) 



with ~ d^/dx^ and = 7*^ <i>,j. We then obtain a 
conservative fluid momentum equation 



dt 



_1 d_ 

y/j dx^ 



V7 



Sfi 



* J dx^ 



(137) 



In the absence of neutrino-matter interactions, and if Carte- 
sian coordinates are used. Equation (137) states that the fluid 
momentum is conserved. 

Using Equations (129), (131), and (134), we can write a 
conservative equation for the energy density Ef + ^ p^ 



- Eii. 



E,+pW + F^ 



where E^, ~ ^ p $ is the gravitational energy density, 



(138) 



F 



d^ 

dx^ 
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(139) 



is the gravitational energy flux density, and $ t = d^/dt. 
Equation (138) states that, in the absence of neutrino-matter 
interactions, the fluid energy (internal plus kinetic) plus the 
gravitational potential energy is conserved. 

7.4. Conservation Laws in Non-Relativistic, Self-Gravitating 
Neutrino Radiation Hydrodynamics 

We discuss conservation laws for non-relativistic, self- 
gravitating neutrino radiation hydrodynamics in this subsec- 
tion. Tracking the evolution of conserved quantities is useful 
when evaluating the physical reliability of numerical simu- 
lations. Equations (123) and (124) (or more precisely, their 
energy-integrated versions) combined with Equations (129)- 
(132) state the conservation of several physical quantities in 
neutrino radiation hydrodynamics. Mass conservation is triv- 
ially stated by Equation (129). Equations (124) and (130) re- 
sult in a total (fluid plus radiation) momentum equation 



-(f. 
dt\ 

1 



Sfl 



-y/7 dx^ 
J dx^ dx^ 



(140) 



Alternatively, combining Equations (124) and (137) results in 



a. . ,.x . 1 s 



dt\ 



H, 



-y/7 dx^ 



sj i 



0. 



(141) 



Equations (123) and (131) result in a conservative equa- 
tion for the "total" (internal plus kinetic plus radiation) energy 
density of a radiating flow 



-pv 



d^ 

dx^ 



(142) 
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Combining Equations (123) and (138) results in a conserva- 
tive equation for the total (internal plus kinetic plus gravita- 
tional plus radiation) energy density of a self-gravitating radi- 
ating flow 



E] 



= 0.(143) 



In Equations (140)-(143), the neutrino energy density, mo- 
mentum density, and stress include contributions from all neu- 
trino species; e.g.. 



(144) 



Equations (141) (assuming Cartesian coordinates are used) 
and (143) are exactly (and locally) conservative. Therefore, 
energy and momentum conservation can be useful checks 
for evaluating the physical reliability of non-relativistic neu- 
trino radiation hydrodynamics simulations based on the 0{v)- 
minus moment equations. However, for radiation hydrody- 
namics based on the C'(u)-plus moment equations, or the 
pseudo-Newtonian moment equations in Section 6, where the 
gravitational potential is obtained by solving a modified Pois- 
son equation (e.g., Kim et al. 2012) and the hydrodynamics 
equations are promoted to include more relativistic effects, we 
are unable to take steps similar to those taken to obtain Equa- 
tions (137) and (138) from Equations (130) and (131), re- 
spectively. Moreover, there are additional source terms on the 
right-hand sides of the total momentum and energy equations 
due to changes in the radiation momentum and energy caused 
by the gravitational field (e.g., the 'bending of light' effect and 
gravitational redshifts). These (locally) non-vanishing grav- 
itational source terms (and the lack of global conservation) 
limit the usefulness of tracking total energy and momentum as 
checks on the physical reliability of (self-gravitating) neutrino 
radiation hydrodynamics simulations. (Note however that in 
3h-1 general relativity, the AMD mass is conserved in asymp- 
totically flat spacetimes; e.g., Baumgarte & Shapiro 2010). 

Consistency with the neutrino number equation results in 
lepton number conservation in (self-gravitating) neutrino ra- 
diation hydrodynamics. Combining the energy-integrated 
version of Equation ( 125) (for electron neutrinos and electron 
antineutrinos) with Equation (132) results in a conservation 
equation for the total electron lepton number 



-f 

dt\ 



Ep, 



F ' 



F 



+ 77,e V 



= 0. (145) 



Lepton number conservation in numerical simulations based 
on solving Equations ( 123) and ( 124), or any of the other radi- 
ation moment equations presented in this paper, may depend 
sensitively on the chosen discretization. In particular, the dis- 
cretized neutrino energy and momentum equations must be 
consistent (in the sense discussed in Sections 7.1 and 7.2) 
with a discretized version of the conservative neutrino number 
equation; Equation (125) in the particular case discussed here. 
When discretizing the monochromatic energy and momentum 
equations, attention should be paid to the cancellations that 
occur when deriving the neutrino number equation from the 
energy and momentum equations, which ideally also occur in 



the discrete limit. To achieve this, the individual terms in the 
moment equations for the neutrino four-momentum should be 
discretized in a coordinated rather than independent fashion. 
Despite the lack of energy conservation due to gravitational 
redshifts discussed above, the total lepton number is still con- 
served. Thus, lepton number conservation may serve as an 
extremely useful gauge on the physical consistency of numer- 
ical simulations of neutrino radiation hydrodynamics (e.g., 
Liebendorfer et al. 2004; Lentz et al. 2012a). 

8. SUMMARY AND DISCUSSION 

In preparation for development of numerical methods for 
multidimensional neutrino radiation hydrodynamics, with 
the eventual goal of simulating the explosion mechanism 
of core-collapse supernovae, we have derived conservative, 
monochromatic general relativistic moment equations for the 
radiation four-momentum (cf. Equation (47)). The radi- 
ation moment equations are conservative in the sense that 
(modulo radiation-matter interactions and geometry sources) 
the time rate of change of the radiation four-momentum is 
governed by space and momentum space divergences. We 
have used the freedom to choose distinct spacetime and 
momentum space coordinates (cf. Lindquist 1966; Riffert 
1986; Mezzacappa & Matzner 1989; Cardall & Mezzacappa 
2003). The evolved radiation quantities are functions of the 
coordinate basis spacetime position components and the 
radiation energy e measured by an observer comoving with 
the fluid. (When integrated over the comoving-frame en- 
ergy, the equations reduce to familiar position space conser- 
vation laws.) The specific choice of phase space coordinates 
is motivated by our intent to develop numerical methods for 
computer simulations of multidimensional neutrino transport. 
Neutrino-matter interactions are on the one hand most eas- 
ily handled computationally in the comoving-frame. On the 
other hand, conservation of global quantities (e.g., energy 
and lepton number) is naturally expressed in the laboratory 
frame. Convenient treatment of neutrino-matter interactions 
and global conservation are naturally handled with the cho- 
sen phase space coordinates. 

We have presented radiation moment equations valid for 
conformally flat spacetimes (Section 5; cf. Equations (67) 
and (72)). We have further specialized the radiation mo- 
ment equations to the pseudo-Newtonian and the Newtonian 
gravity, 0{v) limits (Sections 6 and 7). Furthermore, in the 
0{v) limit, we have presented the ©(wj-plus (Section 7.1) 
and C'(w)-minus (Section 7.2) moment equations. The 0{v)- 
plus radiation energy and momentum equations are given 
by Equations (93) and (104), respectively. In the no grav- 
ity, strict 0{v) limit, these equations are conservative for- 
mulations of similar, non-conservative equations presented 
by other authors (e.g., Buchler 1979; Kaneko et al. 1984; 
Munier & Weaver 1986b). The 0(w)-minus radiation en- 
ergy and momentum equations are given by Equations (123) 
and (124), respectively. Special relativistic moment equations 
are given in Appendix A, while general relativistic moment 
equations for spherically symmetric spacetimes are given in 
Appendix B. The moment equations given in the appen- 
dices are also conservative versions of non-conservative equa- 
tions presented by other authors (e.g.. Castor 1972; Mihalas 
1980; Munier & Weaver 1986b; Muller et al. 2010). (See 
Shibata et al. 201 1; Cardall et al. 2012, for conservative 3h-1 
general relativistic moment equations.) 

We have paid special attention to the issue of neutrino num- 
ber (and lepton number) conservation. For numerical meth- 
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ods based on solving moment equations for the neutrino four- 
momentum, total lepton number conservation will likely serve 
as a very useful gauge on the physical consistency of simula- 
tions of neutrino radiation hydrodynamics. To this end, we 
have exposed the relationship between the equations for the 
neutrino four-momentum and the neutrino number equation 
with multiple examples (cf. Sections 4.3, 5, 6, and 7, and Ap- 
pendix A). The lab-frame radiation number density is related 
to the lab-frame energy and momentum densities by Equa- 
tion (81). The conservative neutrino number equation and 
the equations for the neutrino four-momentum are similarly 
related. The non-relativistic limits of the radiation moment 
equations are not uniquely defined. We obtain consistency 
with the conservative neutrino number equation by adopting 
different orders of v for the energy and momentum equations 
(cf. the 0(ij)-minus and C'(w)-plus approximations in Section 
7). In particular, as was detailed in Section 7.1, when car- 
rying out the steps to obtain the conservative neutrino num- 
ber equation from the conservative energy and momentum 
equations, we observe that terms emanating from the time 
and space derivatives and the geometry sources cancel with 
terms emanating from the energy derivatives. The remaining 
terms constitute the left-hand side of the number equation. 
(The right-hand sides of the equations are similarly related.) 
The energy and momentum equations are therefore consis- 
tent with the conservative number equation. In the 0{v)- 
plus limit, consistency with the conservative 0{v'^) number 
equation is obtained by adopting different orders of v and re- 
taining some higher-order terms in the energy and momen- 
tum equations (accurate to O(u^) and 0{v), respectively). In 
the 0(w)-minus limit, we obtain consistency with the con- 
servative 0{v) number equation by adopting energy and mo- 
mentum equations accurate to 0{v) and 0(1), respectively. 



Ideally, discrete representations of the radiation energy and 
momentum equations can be constructed so that a conserva- 
tive neutrino number equation can be analogously obtained in 
the discrete limit. The discretization is then consistent with 
neutrino number conservation (which is necessary to ensure 
lepton number conservation; Section 7.4). The realization of 
this consistency in a numerical method for neutrino radiation 
transport based on the two-moment model derived here will 
be the focus of a future study. 

General relativistic effects are important in the core- 
collapse supernova environment (e.g., Bmenn et al. 2001; 
Mulleretal. 2012a), and definitive simulations elucidating 
the explosion mechanism of core-collapse supernovae must 
eventually be performed in full general relativity (possibly 
employing multi-energy and multi-angle neutrino transport). 
To this end, the pseudo-Newtonian and the Newtonian, 0{v)- 
plus equations may represent useful self-consistent approxi- 
mations beyond the C'(w)-minus equations. They are closely 
related to the corresponding equations valid for conformally 
flat spacetimes, which again are close to the fully general rel- 
ativistic equations (Shibata et al. 201 1; Cardall et al. 2012). 
It seems plausible that numerical methods developed for the 
self-consistent moment equations presented in this paper — 
which arguably are easier to work with — can be extended in 
steps of increasing degree of complexity to the fully gen- 
eral relativistic case, as has been done in the case of con- 
servative methods for hydrodynamics (e.g.. Font et al. 1994; 
Banyulsetal. 1997). 

This research was supported by the Office of Advanced 
Scientific Computing Research and the Office of Nuclear 
Physics, U.S. Department of Energy. 



APPENDIX 

CONSERVATIVE SPECIAL RELATIVISTIC MOMENT EQUATIONS 

In this appendix we present conservative, multidimensional, monochromatic moment equations valid in special relativity. The 
equations are obtained from the general relativistic moment equations presented in Sections 4 and 5 in the limit of flat spacetime. 
First we list sufficiently general equations to accommodate Cartesian, spherical polar, and cylindrical coordinates. Then we 
specialize the special relativistic moment equations to spherical symmetry and compare with the corresponding equations given 
by Mihalas & Mihalas (1999). Monochromatic, special relativistic moment equations have been presented elsewhere (e.g.. 
Castor 1972; Mihalas 1980; Munier & Weaver 1986b; Mihalas & Mihalas 1999). However, the equations presented in this 
appendix are in fully conservative form. We also present the conservative neutrino number equation, and discuss the relationship 
between the conservative number equation and the conservative radiation energy and momentum equations. 

The invariant line element is 

ds^ = g^^dx^' dx", (Al) 

where, for flat spacetime, the covariant metric tensor is diagonal and given by g^i, = diag[— 1, 1, a^{x^), b^{x^) c^(x^)]. The 
contravariant metric tensor is g'"' = diag[— 1, 1, a~^, b^^ c^^]. The transformation from the orthonormal tetrad basis to the 
global coordinate basis may simply be written as = diag[l, 1, a^^, b^^ c^^], and the corresponding inverse transformation 
is e'jj = diag[l, 1, a, 6 c]; i.e., e"^ e'jj = 6"^. We also have y/—g = = abc. The Lorentz transformation (or boost) from the 
orthonormal comoving basis to the orthonormal (non-comoving) tetrad basis is given by 

where W = ( 1 — )~^/^ is the Lorentz factor In Equation (A2), and vi are three-velocity components defined with 
respect to the orthonormal tetrad basis, and are therefore accented with a bar (vj = = = vi, for i ~ i) — they should not be 
considered spatial components of a four- vector The placement of, and the accent on, the indices of the three-velocity components 
match those of the Lorentz transformation A^^. This results in unambiguous notation when relations between quantities in the 
different reference frames are made explicit (see for example Equations (A16)-(A17) and (A18)-(A20) below). The Lorentz 
transformation from the orthonormal tetrad basis to the orthonormal comoving basis, is obtained by reversing the sign on the 
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three-velocity components; i.e.. 



AO- AO- 
A*- A*- 



W 



» ^ W+1 " 



(A3) 



With the specifications of the transformations between the orthonormal tetrad basis and the coordinate basis, and the Lorentz 
transformations between the orthonormal comoving basis and the orthonormal tetrad basis we can obtain specific expressions 
for the four-velocity of the comoving observer. Equation (58). We can also relate the Eulerian projections of the number-flux 
four-vector, the stress-energy tensor, and the third-order moments in terms of the corresponding Lagrangian projections (Section 
5). In the comoving-frame, the four-velocity of the comoving observer is simply m'' = (1,0 ). The four-velocity of the comoving 
observer in the orthonormal tetrad basis is 

(A4) 



uf" = A'i = W{l,v'), 



while the coordinate basis four-velocity of the comoving observer is 



i'' = e'^.M^"' = (A5) 

where = e\ u'. Note that = u*. Moreover, by contracting Equation (A5) with u^;, we find = (1 — Viv'- . 

The Eulerian projections of the number-flux four-vector can now be expressed in terms of the corresponding Lagrangian 
projections as (cf. Equations (59) and (60)) 



where T-L^ '\s related to the comoving-frame moments by 



- I - w"^ 



(A6) 
(A7) 



(A8) 



The Eulerian projections of the stress-energy tensor are related to the corresponding Lagrangian projections in a similar manner. 
In particular, we have (cf. Equations (61) and (62)) 



E = W'^ J + 2Wv^'H} ^ V, vj JC\ 
r = WW + W^v' J + Vj JC'^ + W v' Vj W , 



(A9) 
(AlO) 
(All) 



In Equations (A9)-(A1 1), the Lagrangian projections W and /C*^ are related to the comoving-frame moments by 



W+1 



W + \ 



v' Vi 6^- + 6\ Vj 



_____ 



(A12) 
(A13) 



Note that JC^^ = fc*^ J', where the rank-two variable Eddington tensor fc*^ is defined in Equation (26). Finally, we express the 
third-order moments in terms of the Lagrangian projections (cf. Equations (65) and (66)) 

- 2'^'= = C^'' + W (v' JC^'' + JC'' + v'^ JC'^ + [v' n'' + v' v'' W + v^- H'^ + v'^ J, (A14) 

where the Lagrangian projections of the third-order tensor is related to the comoving-frame moments by 



+ 



5", + 



W +1 



Vi 6) 5\ + v^ % 5\ 5\ + v" vj^ S\ S\ 



{w + iy 



v' Vi v' Vj 6\ + Vi v^ Vj. S^- + v^ Vj v^ Vf. 5\ 



3 



+ (VF+1)3 ^'"'^'^'^^^'^0^'"''' ^^^^^ 



and = P^*^ J. The rank-three variable Eddington tensor P^'^ is defined in Equation (32). 

Alternatively, we can write the Eulerian projections of the number-flux four- vector in terms of the number- flux four-vector in 
the orthonormal tetrad basis, whose components are expressed in terms of comoving-frame angular moments by 



N^^^wi^Jj^ + ViH^y 
N'^5\n^ + Wv'{jj^ 



w 



(A16) 
(A17) 
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Similarly, The monochromatic radiation energy density, momentum density, and stress in the orthonormal tetrad basis can be 
expressed in terms of the comoving-frame moments as 



w 



w + \ 



w 



l( 



' Vi Vj IC 



(A18) 
(A19) 

(A20) 



Equations (A18) and (A19) correspond to Equations (91.10) and (91.11) in Mihalas & Mihalas (1999) (see also Equations (182) 
and (183) in Munier & Weaver 1986a). Equation (A20) corresponds to Equation (91.12) in Mihalas & Mihalas (1999) (see also 
Equation (184) in Munier & Weaver 1986a). We also have 



- W^'^ = 5% S^- C^'^ + w( 6\ Sl^ + v^ S\ S'^f IC''' + v' Sl (S*^. IC^'' 

^ 3 k \ 3 k 3 k 



{v^ v'' S\ W + v^ v'' S^. W + v' v^ s\n''] + 



^ {S\ 5^ v% + 5\ 5\ v^vj + 5\ v% ] C^'^ 



+W^v'v^ v'' J + 



w ^-\ 



((5^5 v^vj + v^vi (5^-) v'^r^ + (j\ v^'v'^ + v'^vi (5^ v^K}'' + ((5^^ v'^v-^ + v^^ S\ 



'W + 1 

'{i + wy 



v'v^v^ ViW + 



^——^vVv'^ViVjVl^C^K 



(A21) 



We can now list the conservative special relativistic radiation moment equations. 

The conservative, monochromatic special relativistic radiation moment equations are obtained from Equations (67) and (72) 
by setting a = V' = 1 ^nd /3' = 0. The radiation energy equation becomes 



V7 

where the energy space energy flux has been defined as 



(A22) 



= e { ( + 5/ V, +W-^ Qt V, vk -WH t-,) ^ 



5,^ +W- Q^^ Vk - v^ T- 



dt 

\ dv' 1 



5^ +w- Q'^'' Vk 



(A23) 



We have written d In W/dt — Vi dv' / dt and d In W/dx^ = W"^ Vidv'^ / dx^ . The terms proportional to derivatives of three- 
velocity components account for Doppler shifts of the radiation energy spectrum as measured by the comoving observer The 
radiation momentum equation becomes 





dt 1^/7 dx^ 



2" dx 

where we have defined the energy space momentum flux as 



}_d_ 

e2 de 



(A24) 



dv^ 
'dt 



^77 



(A25) 



We specialize the monochromatic radiation energy and momentum equations to spherical polar coordinates; i.e., (x^ , a;^ , ) = 
(r, 9,4>), a = b = r, and c = sin0, and impose spherical symmetry (d/dO = d/d4> = 0). The only nonzero component of the 
fluid three-velocity is v^ ^ vi = v,.- In spherical symmetry, the comoving-frame angular moments are 



(A26) 
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with fi — cos i9. Then, the monochromatic radiation energy and momentum equations become 



9f 
'dt 



J.2 Qj, \ 



]_d_ 



2ttW 



C[/] d^i + Vr j C[f]tld^i) 



and 



dt 



J__9 / 

J.2 Qj. \ 



1 



r'S -- J-IC - 



]_d_ 

£2 de 



e^S'] =2TrW 



' if] - 



C [/] ^ dyU 



respectively. With spherical symmetry imposed, we have 

E = W^ {J + 2vr'H + vlK.) , 
F^W^ {[l + vl]'H + Vr[J + JC]\ 
S = W^ (lC + 2vrH + v^.j) . 



33 



Moreover, we have S.^ = IC^ = /C^^ = ^[J-JC) and 
derivatives can be written explicitly in terms of the comoving-frame moments as 



{All) 

(A28) 

(A29) 
(A30) 
(A31) 

\{J — K.)- The fluxes inside the energy 



{VrH + IC) + I'r ( I'r JC + C) 

vr {J-K.) + {'H-C)] 

J r ) 



dvr 
'dt 



dv,. 
'dt 



Vr in+VrJC) + iK: + VrC) 



Vr {VrH + K.) + {Vr K. + C) 



dVr 

dr 

dVr 

dr 



where we have used 



-Q^^^^W^ (C + 3vrlC + 3vln + v^rJ] 
e 

-^Qil^lQii^lwi[n-c] + vr [J-JC]). 



(A32) 

(A33) 

(A34) 
(A35) 



Equations (A27) and (A28) are monochromatic lab-frame radiation energy and momentum equations in conservative form — valid 
to all orders of Vr- They can be expressed explicitly in terms of the comoving-frame moments, and compared with corresponding 
non-conservative equations in Mihalas & Mihalas (1999); their Equations (95.11) and (95.12). Equation (A27) is equivalent 
to Equation (95.11) plus iv times Equation (95.12) in Mihalas & Mihalas (1999). Similarly, Equation (A28) is equivalent to 
Equation (95. 12) plus Vr times Equation (95. 1 1) in Mihalas & Mihalas (1999). 

The special relativistic moment equations are consistent with the conservative neutrino number equation, which we obtain by 
adding We~^ times Equation (A22) and —We^^ contracted with Equation (A24). An intermediate result is 



dt\e 
1 
e 



e~v'T, 



. py dx^ f , 



Vf| J",; - (S V, - v,^ 



r - v^ s 



1 d 



dv"^ 
'dt 



de 



1 



- v' S', 
ox-i dx'^ 



1 

e Jn 



C [/] dfl. 



(A36) 



(Note that w,, 



-e.) From Equations (A23) and (A25) we have 
1 



- W"^ (e Vi - v^ T,j v, 



dv'^ 
'dt 



Si 



dii 



(A37) 



Thus, by virtue of the cancellation of the terms on the second line with the last term on the left-hand side. Equation (A36) reduces 
to the conservative monochromatic number equation 



dEj^ _l d_ 

dt dx^ 

where the energy space number flux is 



1 d 



1 



C [/] (A38) 
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J-J - 5\ v'' 
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Equation (A38) is a conservative equation for the neutrino number density. In the absence of neutrino-matter interactions, it 
states that the neutrino number is conserved. In carrying out the derivation of the number equation from the radiation energy 
and momentum equations, we observe that the leftover terms emanating from bringing W and Wv"^ inside the space and time 
derivatives in the energy and momentum equations cancel with the leftover terms from bringing e~ ^ inside the energy derivatives. 
The remaining terms constitute the left-hand side of the conservative number equation (Equation ( A38)). These cancellations may 
be key to constructing a numerical scheme for neutrino transport based on the solution of the conservative energy and momentum 
equations that is also consistent with the conservative number equation. Cancellations similar to those occurring in the continuum 
derivation of the number equation form the energy and momentum equations should also occur in the discrete limit in order to 
ensure consistency with neutrino number conservation. Such consistency may help ensure lepton conservation in simulations of 
neutrino radiation hydrodynamics. 

GENERAL RELATIVISTIC MOMENT EQUATIONS FOR SPHERICALLY SYMMETRIC SPACETIMES 

Miiller et al. (2010, hereafter referred to as MID 10) have recently presented numerical methods for multidimensional, general 
relativistic neutrino radiation hydrodynamics for the case where the conformal flatness condition (CFC) is imposed on the spatial 
metric (e.g., Wilson et al. 1996). MJDIO employ the so-called ray-by-ray approach to multidimensional neutrino transport, 
where the radiation flux is assumed to be purely radial, and the radiation field is essentially obtained by solving independent 
spherically symmetric problems along each radial ray. (Lateral advection of neutrinos in optically thick regions and non-radial 
components of the radiation pressure gradient are still taken into account.) The radiation moment equations for the conformally 
flat spacetime become fully general relativistic when spherical symmetry is imposed. For the sake of comparing the conservative 
moment equations presented here with the equations solved by MID 1 0, we adopt a spherically symmetric spacetime, and list — in 
full detail — general relativistic moment equations. The equations presented here can be obtained directly from Equations (67), 
(72), and (83). 

The invariant spacetime line element can be decomposed as 



ds^ = -a'^dt^ + 7y ( dx' + P'dt ) ( dx^ + dt ) . 



(Bl) 



We adopt spherical polar coordinates (x^, x^, x'^) = {r,9,(t)). With spherical symmetry imposed, non-spectral quantities are 
only functions of time t and the radial coordinate r. The only nonzero component of the fluid three-velocity in the orthonormal 
tetrad basis is = vi = v^. For the CFC spacetime, the spatial metric 7^ = is diagonal, and the conformal metric is 

7ij = diag[ 1, , sin^ 6 ], and ?/' is the conformal factor Moreover, a is the lapse function and is the shift vector (only the 
radial component of the shift vector (3^ is nonzero when spherical symmetry is imposed). The lapse function, the shift vector, and 
the conformal factor can be obtained by solving a system of nonlinear elliptic equations (Wilson et al. 1996). For the spherically 
symmetric spacetime, the determinant of the spatial metric is ^/^ = tp^ sin 6, and we can set = [ 1, —/?'', 0, ], 
e%= [0, 0, 0],ande\ = ?A-2(jjjjg[i^ r-\ {r sinO)-^]. 

The monochromatic lab-frame radiation energy and momentum equations, expressed in terms of the comoving-frame angular 
moments J^, H, /C, and C, become 
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(B2) 
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and 
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(B3) 



respectively. Equations (B2) and (B3) are conservative evolution equations for the monochromatic lab-frame radiation energy 
density and momentum density, respectively. The square root of the spatial metric determinant has been absorbed in quantities 
accented with a hat (e.g., J' = y/jJ'', cf. MJDIO). The expression inside the time derivative in Equation (B2) is the contribution 
from radiation to the "matter sources" in the definition of the (conserved) ADM mass (cf. MJDIO). For flat spacetime (i.e., 
a ^ il) = 1 and = 0), Equations (B2) and (B3) reduce to Equations (A27) and (A28), respectively. In Equations (B2) and 
(B3), we have defined the "convective derivative" and the "proper time derivative" along constant coordinate lines; 
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respectively. 

The conservative, monochromatic lab-frame number equation, also expressed in terms of comoving-frame angular moments, 
can be obtained directly from Equation (83), or by adding We^^ times Equation (B2) and —VrWe^^ times Equation (B3). The 
result is 
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where the comoving-frame number density and number flux are 



(B5) 



(B6) 



Note that {J,T-L} = e { Jm, T-Lj^ } (cf. Equation (A26)). Equation (B5) corresponds to the evolution equation for the neutrino 
number density given by MJDIO (their Equation (30); corrected here for a few misprints). For flat spacetime. Equation (B5) 
reduces to the corresponding special relativistic equation (cf. Equation (A38)). 

MJD 10 (cf. their Appendix B) have taken important initial steps in developing numerical methods for neutrino transport based 
on moments models where simultaneous conservation of energy and lepton number is considered (see Liebendorfer et al. 2004, 
for the Boltzmann case in spherical symmetry). There may, however, be room for further improvements. Equations (B2) and 
(B3) differ from the energy and momentum equations solved by MJDIO (their Equations (27) and (28), respectively). We obtain 
their energy equation by adding W times Equation (B2) and —Vr W times Equation (B3). Similarly, we obtain their momentum 
equation by adding —Vr W times Equation (B2) and W times Equation (B3). The relationship between the equation for the 
neutrino number density and the neutrino energy density in MJDIO (their Equations (30) and (27), respectively) is relatively 
simple: they differ only by a factor e. Note that there is a similarly simple relationship between Equation (52) and the conservative 
number equation (Equation (41)). This simple relationship is convenient when constructing a numerical method based on the 
radiation energy and momentum equations, that is simultaneously consistent with the conservative neutrino number equation 
(which ensures number conservation). However, the energy and momentum equations solved by MJDIO (obtainable directly 
from Equation (52)) are formulated in a non-conservative form. The "source terms" depend on time and space derivatives of the 
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fluid three-velocity, which do not vanish in the absence of gravity and neutrino-matter interactions. Conservation of energy is 
not guaranteed when the non-conservative formulation of the energy equation is used. (Additional complications may also arise 
in the presence of shocks.) As a possible improvement, a numerical scheme that simultaneously conserves energy and lepton 
number may be based on solving Equations (B2) and (B3), which are conservative formulations of the radiation energy and 
momentum equations, respectively. However, such a scheme is potentially much more complicated than the algorithm outlined 
in Appendix B in MJDIO. When carrying out the the steps for obtaining the conservative number equation from the conservative 
energy and momentum equations, we observe that terms emanating from the time derivatives, space derivatives, and the geometry 
sources cancel with terms emanating from the energy derivatives. The remaining terms constitute the left-hand side of the number 
equation. Similar cancellations must occur in the discrete limit in order for the discretized energy and momentum to be consistent 
with neutrino number conservation (which is necessary to ensure lepton number conservation). 
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